From 89c63155b96f16e5861e3002a5f54d37b09d119c Mon Sep 17 00:00:00 2001 From: kmtsui Date: Tue, 18 Jul 2023 14:20:00 +0900 Subject: [PATCH] Save ID/OD boundary traj points in WCSimRootTrack --- include/WCSimRootEvent.hh | 20 ++++++++++++--- include/WCSimTrajectory.hh | 25 +++++++++++++++++++ src/WCSimEventAction.cc | 20 ++++++++++++--- src/WCSimRootEvent.cc | 23 +++++++++++++++--- src/WCSimTrajectory.cc | 50 +++++++++++++++++++++++++++++++++++++- 5 files changed, 125 insertions(+), 13 deletions(-) diff --git a/include/WCSimRootEvent.hh b/include/WCSimRootEvent.hh index c86b773ad..44bdc6a8d 100644 --- a/include/WCSimRootEvent.hh +++ b/include/WCSimRootEvent.hh @@ -43,6 +43,9 @@ private: Double_t fTime; Int_t fId; Int_t fParentId; + std::vector> boundaryPoints; + std::vector boundaryKEs; + std::vector boundaryTypes; // 1 = blacksheet, 2 = tyvek, 3 = cave public: WCSimRootTrack() {} @@ -60,7 +63,10 @@ public: Int_t parenttype, Double_t time, Int_t id, - Int_t idParent); + Int_t idParent, + std::vector> bPs, + std::vector bKEs, + std::vector bTypes); virtual ~WCSimRootTrack() { } bool CompareAllVariables(const WCSimRootTrack * c) const; @@ -79,8 +85,11 @@ public: Double_t GetTime() const { return fTime;} Int_t GetId() const {return fId;} Int_t GetParentId() const {return fParentId;} + std::vector> GetBoundaryPoints() {return boundaryPoints;} + std::vector GetBoundaryKEs() {return boundaryKEs;} + std::vector GetBoundaryTypes() {return boundaryTypes;} - ClassDef(WCSimRootTrack,3) + ClassDef(WCSimRootTrack,4) }; @@ -447,7 +456,10 @@ public: Int_t parenttype, Double_t time, Int_t id, - Int_t idParent); + Int_t idParent, + std::vector> bPs, + std::vector bKEs, + std::vector bTypes); WCSimRootTrack * AddTrack (WCSimRootTrack * track); WCSimRootTrack * RemoveTrack(WCSimRootTrack * track); @@ -480,7 +492,7 @@ public: TClonesArray *GetCaptures() const {return fCaptures;} - ClassDef(WCSimRootTrigger,4) //WCSimRootEvent structure + ClassDef(WCSimRootTrigger,5) //WCSimRootEvent structure }; diff --git a/include/WCSimTrajectory.hh b/include/WCSimTrajectory.hh index e9221d69c..9b4213277 100644 --- a/include/WCSimTrajectory.hh +++ b/include/WCSimTrajectory.hh @@ -72,6 +72,26 @@ public: // with description inline void SetStoppingVolume(G4VPhysicalVolume* currentVolume) { stoppingVolume = currentVolume;} +// Functions to Set/Get boundary points + inline void SetBoundaryPoints(std::vector> bPs, + std::vector bKEs, + std::vector bTypes) + { + boundaryPoints = bPs; + boundaryKEs = bKEs; + boundaryTypes = bTypes; + } + inline void AddBoundaryPoint(std::vector bPs, + G4double bKEs, + G4int bTypes) + { + boundaryPoints.push_back(bPs); + boundaryKEs.push_back(bKEs); + boundaryTypes.push_back(bTypes); + } + inline std::vector> GetBoundaryPoints() {return boundaryPoints;} + inline std::vector GetBoundaryKEs() {return boundaryKEs;} + inline std::vector GetBoundaryTypes() {return boundaryTypes;} // Other member functions virtual void ShowTrajectory(std::ostream& os=G4cout) const; @@ -107,6 +127,11 @@ public: // with description G4bool SaveIt; G4String creatorProcess; G4double globalTime; + + // Boundary points; + std::vector> boundaryPoints; + std::vector boundaryKEs; + std::vector boundaryTypes; // 1 = blacksheet, 2 = tyvek, 3 = cave }; /*** TEMP : M FECHNER *********** diff --git a/src/WCSimEventAction.cc b/src/WCSimEventAction.cc index 0b58f949e..2eb9043a8 100644 --- a/src/WCSimEventAction.cc +++ b/src/WCSimEventAction.cc @@ -1210,7 +1210,10 @@ void WCSimEventAction::FillRootEvent(G4int event_id, injhfNtuple.parent[k], injhfNtuple.time[k], 0, - 0); + 0, + std::vector>(), + std::vector(), + std::vector()); } // the rest of the tracks come from WCSimTrajectory @@ -1350,7 +1353,10 @@ void WCSimEventAction::FillRootEvent(G4int event_id, parentType, ttime, id, - idPrnt); + idPrnt, + trj->GetBoundaryPoints(), + trj->GetBoundaryKEs(), + trj->GetBoundaryTypes()); } @@ -1742,7 +1748,10 @@ void WCSimEventAction::FillRootEventHybrid(G4int event_id, injhfNtuple.parent[k], injhfNtuple.time[k], 0, - 0); + 0, + std::vector>(), + std::vector(), + std::vector()); } // the rest of the tracks come from WCSimTrajectory @@ -1878,7 +1887,10 @@ void WCSimEventAction::FillRootEventHybrid(G4int event_id, parentType, ttime, id, - idPrnt); + idPrnt, + trj->GetBoundaryPoints(), + trj->GetBoundaryKEs(), + trj->GetBoundaryTypes()); } diff --git a/src/WCSimRootEvent.cc b/src/WCSimRootEvent.cc index d06e8778c..5dced955a 100644 --- a/src/WCSimRootEvent.cc +++ b/src/WCSimRootEvent.cc @@ -404,7 +404,10 @@ WCSimRootTrack *WCSimRootTrigger::AddTrack(Int_t ipnu, Int_t parenttype, Double_t time, Int_t id, - Int_t idParent) + Int_t idParent, + std::vector> bPs, + std::vector bKEs, + std::vector bTypes) { // Add a new WCSimRootTrack to the list of tracks for this event. // To avoid calling the very time consuming operator new for each track, @@ -428,7 +431,10 @@ WCSimRootTrack *WCSimRootTrigger::AddTrack(Int_t ipnu, parenttype, time, id, - idParent); + idParent, + bPs, + bKEs, + bTypes); fNtrack++; return track; } @@ -466,7 +472,10 @@ WCSimRootTrack *WCSimRootTrigger::AddTrack(WCSimRootTrack * track) track->GetParenttype(), track->GetTime(), track->GetId(), - track->GetParentId()); + track->GetParentId(), + track->GetBoundaryPoints(), + track->GetBoundaryKEs(), + track->GetBoundaryTypes()); fNtrack++; return track_out; } @@ -497,7 +506,10 @@ WCSimRootTrack::WCSimRootTrack(Int_t ipnu, Int_t parenttype, Double_t time, Int_t id, - Int_t idParent) + Int_t idParent, + std::vector> bPs, + std::vector bKEs, + std::vector bTypes) { // Create a WCSimRootTrack object and fill it with stuff @@ -521,6 +533,9 @@ WCSimRootTrack::WCSimRootTrack(Int_t ipnu, fTime = time; fId = id; fParentId = idParent; + boundaryPoints = bPs; + boundaryKEs = bKEs; + boundaryTypes = bTypes; } //_____________________________________________________________________________ diff --git a/src/WCSimTrajectory.cc b/src/WCSimTrajectory.cc index 943887162..23de5996d 100644 --- a/src/WCSimTrajectory.cc +++ b/src/WCSimTrajectory.cc @@ -21,7 +21,11 @@ WCSimTrajectory::WCSimTrajectory() PDGEncoding( 0 ), PDGCharge(0.0), ParticleName(""), initialMomentum( G4ThreeVector() ),SaveIt(false),creatorProcess(""), globalTime(0.0) -{;} +{ + boundaryPoints.clear(); + boundaryKEs.clear(); + boundaryTypes.clear(); +} WCSimTrajectory::WCSimTrajectory(const G4Track* aTrack) { @@ -49,6 +53,10 @@ WCSimTrajectory::WCSimTrajectory(const G4Track* aTrack) } else creatorProcess = ""; + + boundaryPoints.clear(); + boundaryKEs.clear(); + boundaryTypes.clear(); } WCSimTrajectory::WCSimTrajectory(WCSimTrajectory & right):G4VTrajectory() @@ -72,6 +80,10 @@ WCSimTrajectory::WCSimTrajectory(WCSimTrajectory & right):G4VTrajectory() positionRecord->push_back(new G4TrajectoryPoint(*rightPoint)); } globalTime = right.globalTime; + + boundaryPoints = right.boundaryPoints; + boundaryKEs = right.boundaryKEs; + boundaryTypes = right.boundaryTypes; } WCSimTrajectory::~WCSimTrajectory() @@ -84,6 +96,10 @@ WCSimTrajectory::~WCSimTrajectory() positionRecord->clear(); delete positionRecord; + + boundaryPoints.clear(); + boundaryKEs.clear(); + boundaryTypes.clear(); } void WCSimTrajectory::ShowTrajectory(std::ostream& os) const @@ -174,6 +190,32 @@ void WCSimTrajectory::AppendStep(const G4Step* aStep) { positionRecord->push_back( new G4TrajectoryPoint(aStep->GetPostStepPoint()-> GetPosition() )); + + // Save boundary points + G4StepPoint* thePrePoint = aStep->GetPreStepPoint(); + G4VPhysicalVolume* thePrePV = thePrePoint->GetPhysicalVolume(); + + G4StepPoint* thePostPoint = aStep->GetPostStepPoint(); + G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume(); + + if (thePrePV && thePostPV && thePrePV->GetName()!=thePostPV->GetName()) + { + G4String thePrePVName = thePrePV->GetName(); + G4String thePostPVName = thePostPV->GetName(); + G4int ty = 0; + if (thePrePVName.contains("BlackSheet") || thePostPVName.contains("BlackSheet")) ty = 1; + else if (thePrePVName.contains("Cave") || thePostPVName.contains("Cave")) ty = 3; + else if (thePrePVName.contains("Tyvek") || thePostPVName.contains("Tyvek")) ty = 2; + if (ty>0) + { + const G4Track* track = aStep->GetTrack(); + std::vector bPs(3); + bPs[0] = track->GetPosition().x(); bPs[1] = track->GetPosition().y(); bPs[2] = track->GetPosition().z(); + AddBoundaryPoint(bPs, track->GetKineticEnergy(), ty); + // G4cout<<"Step point "<GetCurrentStepNumber () <<" "<GetPosition().x()<<" "<GetPosition().y()<<" "<GetPosition().z()<< + // " "<GetKineticEnergy()<<" "<GetName()<<" "<GetName()<positionRecord)[0]; seco->positionRecord->clear(); + + ent = (seco->GetBoundaryPoints()).size(); + for (G4int i=0;iGetBoundaryPoints()).at(i),(seco->GetBoundaryKEs()).at(i),(seco->GetBoundaryTypes()).at(i)); + } }