10 #include <phparameter/PHParameters.h>
31 #include <calobase/TowerInfo.h>
32 #include <calobase/TowerInfoContainer.h>
33 #include <calobase/TowerInfoContainerv1.h>
34 #include <calobase/TowerInfoDefs.h>
45 #include <Geant4/G4Field.hh>
46 #include <Geant4/G4FieldManager.hh>
47 #include <Geant4/G4ParticleDefinition.hh>
48 #include <Geant4/G4PropagatorInField.hh>
49 #include <Geant4/G4ReferenceCountedHandle.hh>
50 #include <Geant4/G4Step.hh>
51 #include <Geant4/G4StepPoint.hh>
52 #include <Geant4/G4StepStatus.hh>
53 #include <Geant4/G4String.hh>
54 #include <Geant4/G4SystemOfUnits.hh>
55 #include <Geant4/G4ThreeVector.hh>
56 #include <Geant4/G4TouchableHandle.hh>
57 #include <Geant4/G4Track.hh>
58 #include <Geant4/G4TrackStatus.hh>
59 #include <Geant4/G4TransportationManager.hh>
60 #include <Geant4/G4Types.hh>
61 #include <Geant4/G4VPhysicalVolume.hh>
62 #include <Geant4/G4VTouchable.hh>
63 #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"))
114 if (mappingfilename.empty())
121 std::cout <<
PHWHERE <<
" Could not locate " << url << std::endl;
122 std::cout <<
"use empty filename to ignore mapfile" << std::endl;
125 TFile*
file = TFile::Open(url.c_str());
143 catch (std::exception&
e)
145 std::cout << e.what() << std::endl;
156 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
157 G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
182 layer_id = layer_tower.first;
183 tower_id = layer_tower.second;
195 G4StepPoint* prePoint = aStep->GetPreStepPoint();
196 G4StepPoint* postPoint = aStep->GetPostStepPoint();
198 double pretime = prePoint->GetGlobalTime() / nanosecond;
199 double posttime = postPoint->GetGlobalTime() / nanosecond;
200 if (posttime < m_tmin || pretime >
m_tmax)
return false;
201 if ((posttime - pretime) >
m_dt)
return false;
202 G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) /
GeV;
203 const G4Track* aTrack = aStep->GetTrack();
205 double light_yield = eion;
213 const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
214 const G4ThreeVector& worldPosition = postPoint->GetPosition();
215 G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
216 float lx = localPosition.x() /
cm;
217 float lz = fabs(localPosition.z() /
cm);
221 int lcz = (int) (2.0 * lz) + 1;
222 int lcx = (int) (2.0 * (lx + 42.75)) + 1;
237 light_yield = light_yield *
GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
241 unsigned int ieta = tower_id;
242 unsigned int iphi = (
unsigned int) layer_id / 5;
248 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
267 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
268 G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
295 layer_id = layer_tower.first;
296 tower_id = layer_tower.second;
300 layer_id = touch->GetCopyNumber();
303 G4double edep = aStep->GetTotalEnergyDeposit() /
GeV;
304 G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) /
GeV;
305 G4double light_yield = 0;
307 const G4Track* aTrack = aStep->GetTrack();
312 edep = aTrack->GetKineticEnergy() /
GeV;
313 G4Track* killtrack =
const_cast<G4Track*
>(aTrack);
314 killtrack->SetTrackStatus(fStopAndKill);
320 bool geantino =
false;
325 if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
326 aTrack->GetParticleDefinition()->GetParticleName().find(
"geantino") != std::string::npos)
330 G4StepPoint* prePoint = aStep->GetPreStepPoint();
331 G4StepPoint* postPoint = aStep->GetPostStepPoint();
336 switch (prePoint->GetStepStatus())
338 case fPostStepDoItProc:
345 std::cout <<
GetName() <<
": New Hit for " << std::endl;
351 <<
", current trackid: " << aTrack->GetTrackID() << std::endl;
352 std::cout <<
"phys pre vol: " << volume->GetName()
353 <<
" post vol : " << touchpost->GetVolume()->GetName() << std::endl;
354 std::cout <<
" previous phys pre vol: " <<
m_SaveVolPre->GetName()
355 <<
" previous phys post vol: " <<
m_SaveVolPost->GetName() << std::endl;
386 m_Hit->
set_t(0, prePoint->GetGlobalTime() / nanosecond);
405 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
422 std::cout <<
GetName() <<
": hit was not created" << std::endl;
428 <<
", current trackid: " << aTrack->GetTrackID() << std::endl;
429 std::cout <<
"phys pre vol: " << volume->GetName()
430 <<
" post vol : " << touchpost->GetVolume()->GetName() << std::endl;
431 std::cout <<
" previous phys pre vol: " <<
m_SaveVolPre->GetName()
432 <<
" previous phys post vol: " <<
m_SaveVolPost->GetName() << std::endl;
439 std::cout <<
GetName() <<
": hits do not belong to the same track" << std::endl;
441 <<
", current trackid: " << aTrack->GetTrackID()
456 m_Hit->
set_t(1, postPoint->GetGlobalTime() / nanosecond);
471 const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
472 const G4ThreeVector& worldPosition = postPoint->GetPosition();
473 G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
474 float lx = localPosition.x() /
cm;
475 float lz = fabs(localPosition.z() /
cm);
479 int lcz = (int) (2.0 * lz) + 1;
480 int lcx = (int) (2.0 * (lx + 42.75)) + 1;
495 light_yield = light_yield *
GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
507 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
524 if (postPoint->GetStepStatus() == fGeomBoundary ||
525 postPoint->GetStepStatus() == fWorldBoundary ||
526 postPoint->GetStepStatus() == fAtRestDoItProc ||
527 aTrack->GetTrackStatus() == fStopAndKill)
575 m_Hits = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
576 m_AbsorberHits = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
581 std::cout <<
"PHG4OuterHcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
587 std::cout <<
"PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << std::endl;
597 static const std::string h_field_name =
"hOuterHcalField";
601 TH2*
h =
new TH2F(h_field_name.c_str(),
"Magnetic field (Tesla) in HCal;X (cm);Y (cm)", 2400,
602 -300, 300, 2400, -300, 300);
606 std::cout <<
"PHG4OuterHcalSteppingAction::FieldChecker - make a histograme to check outer Hcal field map."
607 <<
" Saved to Fun4AllServer Histo with name " << h_field_name << std::endl;
610 TH2*
h =
dynamic_cast<TH2F*
>(se->
getHisto(h_field_name));
614 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
620 G4ThreeVector globPosition = aStep->GetPreStepPoint()->GetPosition();
622 G4double globPosVec[4] = {0};
623 G4double FieldValueVec[6] = {0};
625 globPosVec[0] = globPosition.x();
626 globPosVec[1] = globPosition.y();
627 globPosVec[2] = globPosition.z();
628 globPosVec[3] = aStep->GetPreStepPoint()->GetGlobalTime();
630 const Int_t binx = h->GetXaxis()->FindBin(globPosVec[0] /
cm);
631 const Int_t biny = h->GetYaxis()->FindBin(globPosVec[1] /
cm);
633 if (h->GetBinContent(binx, binx) == 0)
636 G4TransportationManager* transportMgr =
637 G4TransportationManager::GetTransportationManager();
640 G4PropagatorInField* fFieldPropagator =
641 transportMgr->GetPropagatorInField();
644 G4FieldManager* fieldMgr = fFieldPropagator->FindAndSetFieldManager(volume);
647 const G4Field* pField = fieldMgr->GetDetectorField();
650 pField->GetFieldValue(globPosVec, FieldValueVec);
652 G4ThreeVector FieldValue = G4ThreeVector(FieldValueVec[0],
653 FieldValueVec[1], FieldValueVec[2]);
655 const double B = FieldValue.mag() / tesla;
657 h->SetBinContent(binx, biny, B);
659 std::cout <<
"PHG4OuterHcalSteppingAction::FieldChecker - "
661 <<
", " << biny <<
" := " << B <<
" Tesla @ x,y = " << globPosVec[0] /
cm
662 <<
"," << globPosVec[1] /
cm <<
" cm" << std::endl;
670 nodeItr.
findFirst(
"PHCompositeNode",
"DST"));
673 std::cout <<
"PHComposite node created: DST" << std::endl;