10 #include <phparameter/PHParameters.h>
29 #include <calobase/TowerInfo.h>
30 #include <calobase/TowerInfoContainer.h>
31 #include <calobase/TowerInfoContainerv1.h>
32 #include <calobase/TowerInfoDefs.h>
43 #include <Geant4/G4AffineTransform.hh>
44 #include <Geant4/G4Field.hh>
45 #include <Geant4/G4FieldManager.hh>
46 #include <Geant4/G4LogicalVolume.hh>
47 #include <Geant4/G4NavigationHistory.hh>
48 #include <Geant4/G4ParticleDefinition.hh>
49 #include <Geant4/G4PropagatorInField.hh>
50 #include <Geant4/G4ReferenceCountedHandle.hh>
51 #include <Geant4/G4Step.hh>
52 #include <Geant4/G4StepPoint.hh>
53 #include <Geant4/G4StepStatus.hh>
54 #include <Geant4/G4String.hh>
55 #include <Geant4/G4SystemOfUnits.hh>
56 #include <Geant4/G4ThreeVector.hh>
57 #include <Geant4/G4TouchableHandle.hh>
58 #include <Geant4/G4Track.hh>
59 #include <Geant4/G4TrackStatus.hh>
60 #include <Geant4/G4TransportationManager.hh>
61 #include <Geant4/G4Types.hh>
62 #include <Geant4/G4VPhysicalVolume.hh>
63 #include <Geant4/G4VTouchable.hh>
64 #include <Geant4/G4VUserTrackInformation.hh>
79 , m_Detector(detector)
80 , m_Params(parameters)
81 , m_EnableFieldCheckerFlag(m_Params->get_int_param(
"field_check"))
82 , m_IsActiveFlag(m_Params->get_int_param(
"active"))
83 , m_IsBlackHoleFlag(m_Params->get_int_param(
"blackhole"))
84 , m_NScintiPlates(m_Params->get_int_param(PHG4HcalDefs::
scipertwr) * m_Params->get_int_param(
"n_towers"))
85 , m_LightScintModelFlag(m_Params->get_int_param(
"light_scint_model"))
86 , m_doG4Hit(m_Params->get_int_param(
"saveg4hit"))
87 , m_tmin(m_Params->get_double_param(
"tmin"))
88 , m_tmax(m_Params->get_double_param(
"tmax"))
89 , m_dt(m_Params->get_double_param(
"dt"))
110 for (
int j = 0;
j < 4;
j++)
123 if (mappingfilename.empty())
130 std::cout <<
"use empty filename to ignore mapfile" << std::endl;
134 TFile*
file = TFile::Open(mappingfilename.c_str());
137 for (
int i = 0;
i < 24;
i++)
144 Tilehist +=
"_chimney";
151 std::cout <<
"ERROR: could not find Histogram " << Tilehist << i <<
" in " <<
m_Params->
get_string_param(
"MapFileName") << std::endl;
168 catch (std::exception&
e)
170 std::cout << e.what() << std::endl;
181 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
182 G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
208 sector_id = std::get<0>(layer_tower);
209 layer_id = std::get<1>(layer_tower);
210 tower_id = std::get<2>(layer_tower);
222 G4StepPoint* prePoint = aStep->GetPreStepPoint();
223 G4StepPoint* postPoint = aStep->GetPostStepPoint();
225 double pretime = prePoint->GetGlobalTime() / nanosecond;
226 double posttime = postPoint->GetGlobalTime() / nanosecond;
227 if (posttime < m_tmin || pretime >
m_tmax)
return false;
228 if ((posttime - pretime) >
m_dt)
return false;
229 G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) /
GeV;
230 const G4Track* aTrack = aStep->GetTrack();
232 double light_yield = eion;
240 const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
241 const G4ThreeVector& worldPosition = postPoint->GetPosition();
242 G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
243 float lx = (localPosition.x() /
cm);
244 float ly = (localPosition.y() /
cm);
247 int lcx = (int) (2.0 * lx) + 1;
248 int lcy = (int) (2.0 * (ly + 0.5)) + 1;
250 if ((sector_id == 29) || (sector_id == 30) || (sector_id == 31))
264 if ((lcy >= 1) && (lcy <=
m_MapCorrHist[tower_id]->GetNbinsY()) &&
278 light_yield = light_yield *
GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
282 unsigned int ieta = tower_id;
283 unsigned int iphi = (
unsigned int) layer_id / 5;
289 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
309 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
310 G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
338 sector_id = std::get<0>(layer_tower);
339 layer_id = std::get<1>(layer_tower);
340 tower_id = std::get<2>(layer_tower);
344 layer_id = touch->GetCopyNumber();
348 G4double edep = aStep->GetTotalEnergyDeposit() /
GeV;
349 G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) /
GeV;
350 G4double light_yield = 0;
352 const G4Track* aTrack = aStep->GetTrack();
357 edep = aTrack->GetKineticEnergy() /
GeV;
358 G4Track* killtrack =
const_cast<G4Track*
>(aTrack);
359 killtrack->SetTrackStatus(fStopAndKill);
365 bool geantino =
false;
370 if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
371 aTrack->GetParticleDefinition()->GetParticleName().find(
"geantino") != std::string::npos)
375 G4StepPoint* prePoint = aStep->GetPreStepPoint();
376 G4StepPoint* postPoint = aStep->GetPostStepPoint();
381 switch (prePoint->GetStepStatus())
383 case fPostStepDoItProc:
394 std::cout <<
GetName() <<
": Bad step status combination for the same track " << std::endl;
400 <<
", current trackid: " << aTrack->GetTrackID() << std::endl;
401 std::cout <<
"phys pre vol: " << volume->GetName()
402 <<
" post vol : " << touchpost->GetVolume()->GetName() << std::endl;
403 std::cout <<
" previous phys pre vol: " <<
m_SaveVolPre->GetName()
404 <<
" previous phys post vol: " <<
m_SaveVolPost->GetName() << std::endl;
411 std::cout <<
GetName() <<
": New Hit for " << std::endl;
417 <<
", current trackid: " << aTrack->GetTrackID() << std::endl;
418 std::cout <<
"phys pre vol: " << volume->GetName()
419 <<
" post vol : " << touchpost->GetVolume()->GetName() << std::endl;
420 std::cout <<
" previous phys pre vol: " <<
m_SaveVolPre->GetName()
421 <<
" previous phys post vol: " <<
m_SaveVolPost->GetName() << std::endl;
436 m_Hit->
set_t(0, prePoint->GetGlobalTime() / nanosecond);
456 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
473 std::cout <<
GetName() <<
": hit was not created" << std::endl;
479 <<
", current trackid: " << aTrack->GetTrackID() << std::endl;
480 std::cout <<
"phys pre vol: " << volume->GetName()
481 <<
" post vol : " << touchpost->GetVolume()->GetName() << std::endl;
482 std::cout <<
" previous phys pre vol: " <<
m_SaveVolPre->GetName()
483 <<
" previous phys post vol: " <<
m_SaveVolPost->GetName() << std::endl;
490 std::cout <<
GetName() <<
": hits do not belong to the same track" << std::endl;
492 <<
", current trackid: " << aTrack->GetTrackID()
507 m_Hit->
set_t(1, postPoint->GetGlobalTime() / nanosecond);
522 const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
523 const G4ThreeVector& worldPosition = postPoint->GetPosition();
524 G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
525 float lx = (localPosition.x() /
cm);
526 float ly = (localPosition.y() /
cm);
529 int lcx = (int) (2.0 * lx) + 1;
530 int lcy = (int) (2.0 * (ly + 0.5)) + 1;
532 if ((sector_id == 29) || (sector_id == 30) || (sector_id == 31))
546 if ((lcy >= 1) && (lcy <=
m_MapCorrHist[tower_id]->GetNbinsY()) &&
560 light_yield = light_yield *
GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
573 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
590 if (postPoint->GetStepStatus() == fGeomBoundary ||
591 postPoint->GetStepStatus() == fWorldBoundary ||
592 postPoint->GetStepStatus() == fAtRestDoItProc ||
593 aTrack->GetTrackStatus() == fStopAndKill)
633 std::cout <<
"PHG4OHCalSteppingAction::SetTopNode - unable to find " <<
m_HitNodeName << std::endl;
635 if (!m_AbsorberHitContainer)
639 std::cout <<
"PHG4OHcalSteppingAction::SetTopNode - unable to find " <<
m_AbsorberNodeName << std::endl;
649 static const std::string h_field_name =
"hOHCalField";
653 TH2F*
h =
new TH2F(h_field_name.c_str(),
"Magnetic field (Tesla) in HCal;X (cm);Y (cm)", 2400,
654 -300, 300, 2400, -300, 300);
658 std::cout <<
"PHG4OHCalSteppingAction::FieldChecker - make a histograme to check outer Hcal field map."
659 <<
" Saved to Fun4AllServer Histo with name " << h_field_name << std::endl;
662 TH2F*
h =
dynamic_cast<TH2F*
>(se->
getHisto(h_field_name));
666 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
672 G4ThreeVector globPosition = aStep->GetPreStepPoint()->GetPosition();
674 G4double globPosVec[4] = {0};
675 G4double FieldValueVec[6] = {0};
677 globPosVec[0] = globPosition.x();
678 globPosVec[1] = globPosition.y();
679 globPosVec[2] = globPosition.z();
680 globPosVec[3] = aStep->GetPreStepPoint()->GetGlobalTime();
682 const Int_t binx = h->GetXaxis()->FindBin(globPosVec[0] /
cm);
683 const Int_t biny = h->GetYaxis()->FindBin(globPosVec[1] /
cm);
685 if (h->GetBinContent(binx, binx) == 0)
688 G4TransportationManager* transportMgr = G4TransportationManager::GetTransportationManager();
691 G4PropagatorInField* fFieldPropagator = transportMgr->GetPropagatorInField();
694 G4FieldManager* fieldMgr = fFieldPropagator->FindAndSetFieldManager(volume);
697 const G4Field* pField = fieldMgr->GetDetectorField();
700 pField->GetFieldValue(globPosVec, FieldValueVec);
702 G4ThreeVector FieldValue = G4ThreeVector(FieldValueVec[0],
703 FieldValueVec[1], FieldValueVec[2]);
705 const double Bz = FieldValue.z() / tesla;
707 h->SetBinContent(binx, biny, Bz);
709 std::cout <<
"PHG4OHCalSteppingAction::FieldChecker - "
710 <<
"volume " << volume->GetName() <<
" / " << volume->GetLogicalVolume()->GetName()
712 <<
", " << biny <<
" : Bz= " << Bz <<
" B = " << FieldValue.mag() / tesla
713 <<
" Tesla @ x,y = " << globPosVec[0] /
cm
714 <<
"," << globPosVec[1] /
cm <<
" cm" << std::endl;
725 else if (type ==
"G4HIT_ABSORBER")
730 std::cout <<
"Invalid output hit node type " << type << std::endl;
739 nodeItr.
findFirst(
"PHCompositeNode",
"DST"));
742 std::cout <<
"PHComposite node created: DST" << std::endl;