File indexing completed on 2025-02-23 09:22:43
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032 #include "RE01Trajectory.hh"
0033
0034 #include "RE01TrackInformation.hh"
0035
0036 #include "G4AttDef.hh"
0037 #include "G4AttDefStore.hh"
0038 #include "G4AttValue.hh"
0039 #include "G4Circle.hh"
0040 #include "G4Colour.hh"
0041 #include "G4DynamicParticle.hh"
0042 #include "G4ParticleTable.hh"
0043 #include "G4ParticleTypes.hh"
0044 #include "G4Polyline.hh"
0045 #include "G4PrimaryParticle.hh"
0046 #include "G4UIcommand.hh"
0047 #include "G4UnitsTable.hh"
0048 #include "G4VVisManager.hh"
0049 #include "G4VisAttributes.hh"
0050
0051 G4ThreadLocal G4Allocator<RE01Trajectory>* myTrajectoryAllocator = 0;
0052
0053
0054 RE01Trajectory::RE01Trajectory(const G4Track* aTrack)
0055 : G4VTrajectory(), fPositionRecord(0), fParticleDefinition(0)
0056 {
0057 fParticleDefinition = aTrack->GetDefinition();
0058 fParticleName = fParticleDefinition->GetParticleName();
0059 fPDGCharge = fParticleDefinition->GetPDGCharge();
0060 fPDGEncoding = fParticleDefinition->GetPDGEncoding();
0061 if (fParticleName == "unknown") {
0062 G4PrimaryParticle* pp = aTrack->GetDynamicParticle()->GetPrimaryParticle();
0063 if (pp) {
0064 if (pp->GetCharge() < DBL_MAX) fPDGCharge = pp->GetCharge();
0065 fPDGEncoding = pp->GetPDGcode();
0066 if (pp->GetG4code() != 0) {
0067 fParticleName += " : ";
0068 fParticleName += pp->GetG4code()->GetParticleName();
0069 }
0070 }
0071 }
0072 fTrackID = aTrack->GetTrackID();
0073 RE01TrackInformation* trackInfo = (RE01TrackInformation*)(aTrack->GetUserInformation());
0074 fTrackStatus = trackInfo->GetTrackingStatus();
0075 if (fTrackStatus == 1) {
0076 fParentID = aTrack->GetParentID();
0077 }
0078 else if (fTrackStatus == 2) {
0079 fParentID = trackInfo->GetSourceTrackID();
0080 }
0081 else {
0082 fParentID = -1;
0083 }
0084 fPositionRecord = new RE01TrajectoryPointContainer();
0085 fPositionRecord->push_back(new G4TrajectoryPoint(aTrack->GetPosition()));
0086 fMomentum = aTrack->GetMomentum();
0087 fVertexPosition = aTrack->GetPosition();
0088 fGlobalTime = aTrack->GetGlobalTime();
0089 }
0090
0091
0092 RE01Trajectory::~RE01Trajectory()
0093 {
0094 size_t i;
0095 for (i = 0; i < fPositionRecord->size(); i++) {
0096 delete (*fPositionRecord)[i];
0097 }
0098 fPositionRecord->clear();
0099
0100 delete fPositionRecord;
0101 }
0102
0103
0104 void RE01Trajectory::ShowTrajectory(std::ostream& os) const
0105 {
0106 os << G4endl << "TrackID =" << fTrackID << " : ParentID=" << fParentID
0107 << " : TrackStatus=" << fTrackStatus << G4endl;
0108 os << "Particle name : " << fParticleName << " PDG code : " << fPDGEncoding
0109 << " Charge : " << fPDGCharge << G4endl;
0110 os << "Original momentum : " << G4BestUnit(fMomentum, "Energy") << G4endl;
0111 os << "Vertex : " << G4BestUnit(fVertexPosition, "Length")
0112 << " Global time : " << G4BestUnit(fGlobalTime, "Time") << G4endl;
0113 os << " Current trajectory has " << fPositionRecord->size() << " points." << G4endl;
0114
0115 for (size_t i = 0; i < fPositionRecord->size(); i++) {
0116 G4TrajectoryPoint* aTrajectoryPoint = (G4TrajectoryPoint*)((*fPositionRecord)[i]);
0117 os << "Point[" << i << "]"
0118 << " Position= " << aTrajectoryPoint->GetPosition() << G4endl;
0119 }
0120 }
0121
0122
0123 void RE01Trajectory::DrawTrajectory() const
0124 {
0125 G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
0126 G4ThreeVector pos;
0127
0128 G4Polyline pPolyline;
0129 for (size_t i = 0; i < fPositionRecord->size(); i++) {
0130 G4TrajectoryPoint* aTrajectoryPoint = (G4TrajectoryPoint*)((*fPositionRecord)[i]);
0131 pos = aTrajectoryPoint->GetPosition();
0132 pPolyline.push_back(pos);
0133 }
0134
0135 G4Colour colour(0.2, 0.2, 0.2);
0136 if (fParticleDefinition == G4Gamma::GammaDefinition())
0137 colour = G4Colour(0., 0., 1.);
0138 else if (fParticleDefinition == G4Electron::ElectronDefinition()
0139 || fParticleDefinition == G4Positron::PositronDefinition())
0140 colour = G4Colour(1., 1., 0.);
0141 else if (fParticleDefinition == G4MuonMinus::MuonMinusDefinition()
0142 || fParticleDefinition == G4MuonPlus::MuonPlusDefinition())
0143 colour = G4Colour(0., 1., 0.);
0144 else if (fParticleDefinition->GetParticleType() == "meson") {
0145 if (fPDGCharge != 0.)
0146 colour = G4Colour(1., 0., 0.);
0147 else
0148 colour = G4Colour(0.5, 0., 0.);
0149 }
0150 else if (fParticleDefinition->GetParticleType() == "baryon") {
0151 if (fPDGCharge != 0.)
0152 colour = G4Colour(0., 1., 1.);
0153 else
0154 colour = G4Colour(0., 0.5, 0.5);
0155 }
0156
0157 G4VisAttributes attribs(colour);
0158 pPolyline.SetVisAttributes(attribs);
0159 if (pVVisManager) pVVisManager->Draw(pPolyline);
0160 }
0161
0162
0163 const std::map<G4String, G4AttDef>* RE01Trajectory::GetAttDefs() const
0164 {
0165 G4bool isNew;
0166 std::map<G4String, G4AttDef>* store = G4AttDefStore::GetInstance("RE01Trajectory", isNew);
0167 if (isNew) {
0168 G4String id("ID");
0169 (*store)[id] = G4AttDef(id, "Track ID", "Bookkeeping", "", "G4int");
0170
0171 G4String pid("PID");
0172 (*store)[pid] = G4AttDef(pid, "Parent ID", "Bookkeeping", "", "G4int");
0173
0174 G4String status("Status");
0175 (*store)[status] = G4AttDef(status, "Track Status", "Bookkeeping", "", "G4int");
0176
0177 G4String pn("PN");
0178 (*store)[pn] = G4AttDef(pn, "Particle Name", "Bookkeeping", "", "G4String");
0179
0180 G4String ch("Ch");
0181 (*store)[ch] = G4AttDef(ch, "Charge", "Physics", "e+", "G4double");
0182
0183 G4String pdg("PDG");
0184 (*store)[pdg] = G4AttDef(pdg, "PDG Encoding", "Bookkeeping", "", "G4int");
0185
0186 G4String imom("IMom");
0187 (*store)[imom] = G4AttDef(imom, "Momentum of track at start of trajectory", "Physics",
0188 "G4BestUnit", "G4ThreeVector");
0189
0190 G4String imag("IMag");
0191 (*store)[imag] = G4AttDef(imag, "Magnitude of momentum of track at start of trajectory",
0192 "Physics", "G4BestUnit", "G4double");
0193
0194 G4String vtxPos("VtxPos");
0195 (*store)[vtxPos] =
0196 G4AttDef(vtxPos, "Vertex position", "Physics", "G4BestUnit", "G4ThreeVector");
0197
0198 G4String ntp("NTP");
0199 (*store)[ntp] = G4AttDef(ntp, "No. of points", "Bookkeeping", "", "G4int");
0200 }
0201 return store;
0202 }
0203
0204
0205 std::vector<G4AttValue>* RE01Trajectory::CreateAttValues() const
0206 {
0207 std::vector<G4AttValue>* values = new std::vector<G4AttValue>;
0208
0209 values->push_back(G4AttValue("ID", G4UIcommand::ConvertToString(fTrackID), ""));
0210
0211 values->push_back(G4AttValue("PID", G4UIcommand::ConvertToString(fParentID), ""));
0212
0213 values->push_back(G4AttValue("Status", G4UIcommand::ConvertToString(fTrackStatus), ""));
0214
0215 values->push_back(G4AttValue("PN", fParticleName, ""));
0216
0217 values->push_back(G4AttValue("Ch", G4UIcommand::ConvertToString(fPDGCharge), ""));
0218
0219 values->push_back(G4AttValue("PDG", G4UIcommand::ConvertToString(fPDGEncoding), ""));
0220
0221 values->push_back(G4AttValue("IMom", G4BestUnit(fMomentum, "Energy"), ""));
0222
0223 values->push_back(G4AttValue("IMag", G4BestUnit(fMomentum.mag(), "Energy"), ""));
0224
0225 values->push_back(G4AttValue("VtxPos", G4BestUnit(fVertexPosition, "Length"), ""));
0226
0227 values->push_back(G4AttValue("NTP", G4UIcommand::ConvertToString(GetPointEntries()), ""));
0228
0229 return values;
0230 }
0231
0232
0233 void RE01Trajectory::AppendStep(const G4Step* aStep)
0234 {
0235 fPositionRecord->push_back(new G4TrajectoryPoint(aStep->GetPostStepPoint()->GetPosition()));
0236 }
0237
0238
0239 void RE01Trajectory::MergeTrajectory(G4VTrajectory* secondTrajectory)
0240 {
0241 if (!secondTrajectory) return;
0242
0243 RE01Trajectory* seco = (RE01Trajectory*)secondTrajectory;
0244 G4int ent = seco->GetPointEntries();
0245
0246
0247 for (int i = 1; i < ent; i++) {
0248 fPositionRecord->push_back((*(seco->fPositionRecord))[i]);
0249 }
0250 delete (*seco->fPositionRecord)[0];
0251 seco->fPositionRecord->clear();
0252 }
