1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | #include "GlueXSteppingAction.hh" |
8 | #include "GlueXUserOptions.hh" |
9 | #include "GlueXPrimaryGeneratorAction.hh" |
10 | #include "GlueXUserEventInformation.hh" |
11 | #include "GlueXUserTrackInformation.hh" |
12 | #include "G4ParallelWorldProcess.hh" |
13 | #include "G4SystemOfUnits.hh" |
14 | #include "G4RunManager.hh" |
15 | #include "G4Threading.hh" |
16 | #include "G4Event.hh" |
17 | #include "G4Decay.hh" |
18 | |
19 | #include <stdio.h> |
20 | #include <exception> |
21 | #include <map> |
22 | |
23 | |
24 | #if BACKGROUND_PROFILING |
25 | #include <TFile.h> |
26 | #include <TTree.h> |
27 | TFile *bgprofiles_file = 0; |
28 | std::map<int, TTree*> bgprofiles; |
29 | #endif |
30 | |
31 | G4Mutex GlueXSteppingAction::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
32 | int GlueXSteppingAction::fStopTracksInCollimator = 0; |
33 | int GlueXSteppingAction::fSaveTrajectories = 0; |
34 | |
35 | GlueXSteppingAction::GlueXSteppingAction() |
36 | { |
37 | G4AutoLock barrier(&fMutex); |
38 | GlueXUserOptions *user_opts = GlueXUserOptions::GetInstance(); |
39 | if (user_opts == 0) { |
40 | G4cerr(*G4cerr_p) << "Error in GlueXSteppingAction constructor - " |
41 | << "GlueXUserOptions::GetInstance() returns null, " |
42 | << "cannot continue." << G4endlstd::endl; |
43 | exit(-1); |
44 | } |
45 | |
46 | std::map<int,double> showersInCol; |
47 | if (user_opts->Find("SHOWERSINCOL", showersInCol)) { |
48 | fStopTracksInCollimator = (showersInCol[1] == 0); |
49 | } |
50 | else { |
51 | fStopTracksInCollimator = 0; |
52 | } |
53 | |
54 | std::map<int,int> trajectories; |
55 | if (user_opts->Find("TRAJECTORIES", trajectories)) { |
56 | fSaveTrajectories = trajectories[1]; |
57 | } |
58 | else { |
59 | fSaveTrajectories = 0; |
60 | } |
61 | |
62 | #if BACKGROUND_PROFILING |
63 | if (bgprofiles_file == 0) { |
64 | bgprofiles_file = new TFile("bgprofiles.root", "recreate"); |
65 | } |
66 | #endif |
67 | |
68 | } |
69 | |
70 | GlueXSteppingAction::GlueXSteppingAction(const GlueXSteppingAction &src) |
71 | {} |
72 | |
73 | GlueXSteppingAction::~GlueXSteppingAction() |
74 | { |
75 | #if BACKGROUND_PROFILING |
76 | G4AutoLock barrier(&fMutex); |
77 | if (bgprofiles_file) { |
78 | std::map<int, TTree*>::iterator iter; |
79 | for (iter = bgprofiles.begin(); iter != bgprofiles.end(); ++iter) { |
80 | iter->second->Write(); |
81 | } |
82 | delete bgprofiles_file; |
83 | bgprofiles_file = 0; |
84 | } |
85 | #endif |
86 | } |
87 | |
88 | void GlueXSteppingAction::UserSteppingAction(const G4Step* step) |
89 | { |
90 | G4Track *track = (G4Track*)step->GetTrack(); |
91 | const G4Event *event = G4RunManager::GetRunManager() |
92 | ->GetCurrentEvent(); |
93 | GlueXUserEventInformation *eventinfo = (GlueXUserEventInformation*) |
| 3 | | 'eventinfo' initialized here | |
|
94 | event->GetUserInformation(); |
| 1 | Calling 'G4Event::GetUserInformation' | |
|
| 2 | | Returning from 'G4Event::GetUserInformation' | |
|
95 | GlueXUserTrackInformation *trackinfo; |
96 | trackinfo = (GlueXUserTrackInformation*)track->GetUserInformation(); |
97 | const G4Step *hyperstep = G4ParallelWorldProcess::GetHyperStep(); |
98 | G4VPhysicalVolume *pvol = hyperstep->GetPostStepPoint()->GetPhysicalVolume(); |
99 | |
100 | |
101 | |
102 | if (fStopTracksInCollimator) { |
| 4 | | Assuming 'fStopTracksInCollimator' is 0 | |
|
| |
103 | if (pvol && (pvol->GetName() == "INSU" || |
104 | pvol->GetName() == "PCTT" || |
105 | pvol->GetName() == "PCPB") ) |
106 | { |
107 | track->SetTrackStatus(fStopAndKill); |
108 | } |
109 | } |
110 | |
111 | |
112 | |
113 | if (pvol && (pvol->GetName() == "World")) { |
114 | track->SetTrackStatus(fStopAndKill); |
115 | } |
116 | |
117 | |
118 | if (trackinfo) { |
| 6 | | Assuming 'trackinfo' is non-null | |
|
| |
119 | int primeID = trackinfo->GetGlueXTrackID(); |
120 | if (primeID > 0) { |
| 8 | | Assuming 'primeID' is > 0 | |
|
| |
121 | const G4VProcess* process; |
122 | process = step->GetPostStepPoint()->GetProcessDefinedStep(); |
123 | if (process && process->GetProcessType() == fDecay) { |
| |
124 | G4TrackVector &secondary = *(G4TrackVector*) |
125 | step->GetSecondaryInCurrentStep(); |
126 | G4TrackVector::iterator iter; |
127 | for (iter = secondary.begin(); iter != secondary.end(); ++iter) { |
| 11 | | Loop condition is false. Execution continues on line 135 | |
|
128 | int newID = eventinfo->AssignNextGlueXTrackID(); |
129 | trackinfo = new GlueXUserTrackInformation(); |
130 | if (eventinfo) { |
131 | trackinfo->SetGlueXTrackID(newID); |
132 | } |
133 | (*iter)->SetUserInformation(trackinfo); |
134 | } |
135 | if (eventinfo) { |
| 12 | | Assuming 'eventinfo' is null | |
|
| |
136 | int mech[2]; |
137 | char *cmech = (char*)mech; |
138 | snprintf(cmech, 5, "%c%c%c%c", 'D', 'C', 'A', 'Y'); |
139 | eventinfo->AddSecondaryVertex(secondary, primeID, mech[0]); |
140 | } |
141 | } |
142 | } |
143 | } |
144 | |
145 | |
146 | if (fSaveTrajectories) { |
| 14 | | Assuming 'fSaveTrajectories' is not equal to 0 | |
|
| |
147 | eventinfo->AddMCtrajectoryPoint(*step, fSaveTrajectories); |
| 16 | | Called C++ object pointer is null |
|
148 | } |
149 | |
150 | #if BACKGROUND_PROFILING |
151 | |
152 | |
153 | |
154 | |
155 | |
156 | |
157 | |
158 | |
159 | |
160 | struct profiler_row_t { |
161 | float totE; |
162 | float x[7]; |
163 | float ppol; |
164 | float xspot[2]; |
165 | int ptype; |
166 | int det; |
167 | int mint; |
168 | float xint[9][3]; |
169 | }; |
170 | const int mint_max = 9; |
171 | static struct profiler_row_t prow[256]; |
172 | int id = G4Threading::G4GetThreadId() + 1; |
173 | assert(id < 256)((id < 256) ? static_cast<void> (0) : __assert_fail ( "id < 256", "src/GlueXSteppingAction.cc", 173, __PRETTY_FUNCTION__ )); |
174 | |
175 | if (track->GetCurrentStepNumber() == 1) { |
176 | if (track->GetParentID() == 0) { |
177 | G4StepPoint *point = step->GetPostStepPoint(); |
178 | const G4ThreeVector &pos = point->GetPosition(); |
179 | prow[id].xspot[0] = pos[0]/cm; |
180 | prow[id].xspot[1] = pos[1]/cm; |
181 | prow[id].mint = 0; |
182 | } |
183 | else { |
184 | G4StepPoint *point = step->GetPreStepPoint(); |
185 | const G4ThreeVector &pos = point->GetPosition(); |
186 | prow[id].xint[prow[id].mint][0] = pos[0]/cm; |
187 | prow[id].xint[prow[id].mint][1] = pos[1]/cm; |
188 | prow[id].xint[prow[id].mint][2] = pos[2]/cm; |
189 | if (prow[id].mint < mint_max - 1) |
190 | ++prow[id].mint; |
191 | } |
192 | } |
193 | |
194 | int det=0; |
195 | if (pvol == 0) |
196 | det = 0; |
197 | else if (pvol->GetName() == "DET1") |
198 | det = 1; |
199 | else if (pvol->GetName() == "DET2") |
200 | det = 2; |
201 | else if (pvol->GetName() == "DET3") |
202 | det = 3; |
203 | else if (pvol->GetName() == "DET4") |
204 | det = 4; |
205 | else if (pvol->GetName() == "DET5") |
206 | det = 5; |
207 | else if (pvol->GetName() == "DET6") |
208 | det = 6; |
209 | else if (pvol->GetName() == "DET7") |
210 | det = 7; |
211 | else if (pvol->GetName() == "DET8") |
212 | det = 8; |
213 | else |
214 | det = 0; |
215 | |
216 | TTree *proftree; |
217 | if (det > 0) { |
218 | G4AutoLock barrier(&fMutex); |
219 | try { |
220 | proftree = bgprofiles.at(det); |
221 | } |
222 | catch(std::out_of_range &err) { |
223 | std::stringstream names; |
224 | names << "det" << det; |
225 | std::stringstream titles; |
226 | titles << "hits in virtual detector " << det; |
227 | proftree = new TTree(names.str().c_str(), titles.str().c_str()); |
228 | proftree->Branch("totE", &prow[0].totE, "totE/F"); |
229 | proftree->Branch("x", &prow[0].x[0], "x[7]/F"); |
230 | proftree->Branch("ppol", &prow[0].ppol, "ppol/F"); |
231 | proftree->Branch("xspot", &prow[0].xspot[0], "xspot[2]/F"); |
232 | proftree->Branch("ptype", &prow[0].ptype, "ptype/I"); |
233 | proftree->Branch("det", &prow[0].det, "det/I"); |
234 | proftree->Branch("mint", &prow[0].mint, "mint/I"); |
235 | proftree->Branch("xint", &prow[0].xint[0][0], "xint[mint][3]/F"); |
236 | bgprofiles[det] = proftree; |
237 | } |
238 | G4StepPoint *point = step->GetPostStepPoint(); |
239 | const G4ThreeVector &pos = point->GetPosition(); |
240 | const G4ThreeVector &mom = point->GetMomentum(); |
241 | const G4ThreeVector &pol = point->GetPolarization(); |
242 | double pmag = mom.mag(); |
243 | double Etot = point->GetTotalEnergy(); |
244 | int pdgcode = track->GetDynamicParticle()->GetPDGcode(); |
245 | int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgcode); |
246 | prow[0].totE = Etot/GeV; |
247 | prow[0].x[0] = pos[0]/cm; |
248 | prow[0].x[1] = pos[1]/cm; |
249 | prow[0].x[2] = pos[2]/cm; |
250 | prow[0].x[3] = mom[0]/pmag; |
251 | prow[0].x[4] = mom[1]/pmag; |
252 | prow[0].x[5] = mom[2]/pmag; |
253 | prow[0].x[6] = pmag/GeV; |
254 | prow[0].ppol= pol.mag(); |
255 | prow[0].xspot[0] = prow[id].xspot[0]; |
256 | prow[0].xspot[1] = prow[id].xspot[1]; |
257 | prow[0].ptype = g3type; |
258 | prow[0].det = det; |
259 | prow[0].mint = prow[id].mint; |
260 | for (int i=0; i < prow[0].mint; ++i) { |
261 | prow[0].xint[i][0] = prow[id].xint[i][0]; |
262 | prow[0].xint[i][1] = prow[id].xint[i][1]; |
263 | prow[0].xint[i][2] = prow[id].xint[i][2]; |
264 | } |
265 | proftree->Fill(); |
266 | } |
267 | |
268 | #endif |
269 | |
270 | } |