1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: OlapEventAction.cc,v 1.5 2006/06/29 17:22:54 gunter Exp $ |
---|
28 | // GEANT4 tag $Name: $ |
---|
29 | // |
---|
30 | // |
---|
31 | // -------------------------------------------------------------- |
---|
32 | // OlapEventAction |
---|
33 | // |
---|
34 | // Author: Martin Liendl - Martin.Liendl@cern.ch |
---|
35 | // |
---|
36 | // -------------------------------------------------------------- |
---|
37 | // |
---|
38 | #include "globals.hh" |
---|
39 | #include <vector> |
---|
40 | #include <sstream> |
---|
41 | |
---|
42 | #include "G4VVisManager.hh" |
---|
43 | #include "G4Trajectory.hh" |
---|
44 | #include "G4TrajectoryContainer.hh" |
---|
45 | #include "G4UImanager.hh" |
---|
46 | #include "G4Event.hh" |
---|
47 | #include "G4RunManager.hh" |
---|
48 | #include "G4Run.hh" |
---|
49 | |
---|
50 | #include "OlapManager.hh" |
---|
51 | #include "OlapEventAction.hh" |
---|
52 | #include "OlapRunAction.hh" |
---|
53 | #include "OlapDetConstr.hh" |
---|
54 | #include "OlapGenerator.hh" |
---|
55 | |
---|
56 | //#define OLAP_DEBUG |
---|
57 | //#define OLAP_DEBUG_2 |
---|
58 | |
---|
59 | std::ostream & |
---|
60 | operator<<(std::ostream& flux, OlapInfo & oi) |
---|
61 | { |
---|
62 | OlapStepInfo si; |
---|
63 | si.thePoint = oi.v1; |
---|
64 | si.theHist = oi.hist1; |
---|
65 | flux << "vol 1: point=" << oi.v1 << " vol 2: point=" << oi.v2 |
---|
66 | << " axis=" << oi.axis << " [" << oi.probNot << "] " |
---|
67 | << oi.info << G4endl; |
---|
68 | if (oi.stAB.size()) |
---|
69 | { |
---|
70 | flux << oi.stAB << G4endl << G4endl; |
---|
71 | flux << oi.stBA << G4endl; |
---|
72 | } |
---|
73 | flux << "A -> B:" << G4endl; |
---|
74 | flux << si << G4endl; |
---|
75 | flux << "B -> A:" << G4endl; |
---|
76 | si.thePoint=oi.v2; |
---|
77 | si.theHist=oi.hist2; |
---|
78 | flux << si << G4endl; |
---|
79 | return flux; |
---|
80 | } |
---|
81 | |
---|
82 | |
---|
83 | OlapInfo::~OlapInfo() |
---|
84 | { |
---|
85 | while( !stAB.empty() ) |
---|
86 | { |
---|
87 | delete stAB.back(); |
---|
88 | stAB.pop_back(); |
---|
89 | } |
---|
90 | |
---|
91 | while( !stBA.empty() ) |
---|
92 | { |
---|
93 | delete stBA.back(); |
---|
94 | stBA.pop_back(); |
---|
95 | } |
---|
96 | |
---|
97 | } |
---|
98 | |
---|
99 | G4bool OlapInfo::operator==(const OlapInfo & rh) |
---|
100 | { |
---|
101 | /* |
---|
102 | const G4AffineTransform & in1 = rh.hist1.GetTopTransform(); |
---|
103 | const G4AffineTransform & in2 = rh.hist2.GetTopTransform(); |
---|
104 | const G4AffineTransform & th1 = hist1.GetTopTransform(); |
---|
105 | const G4AffineTransform & th2 = hist2.GetTopTransform(); |
---|
106 | |
---|
107 | G4cout << "histories: 1 2: " << G4endl; |
---|
108 | for (G4int i=0; i<16; i++) { |
---|
109 | G4cout << th1[i] << '\t' << th2[i] << G4endl; |
---|
110 | G4cout << in1[i] << '\t' << in2[i] << G4endl; |
---|
111 | G4cout << "----------------------------" << G4endl; |
---|
112 | } |
---|
113 | */ |
---|
114 | G4bool result = false; |
---|
115 | if ( |
---|
116 | ( |
---|
117 | ( hist1.GetTopVolume() == rh.hist1.GetTopVolume() ) && |
---|
118 | ( hist2.GetTopVolume() == rh.hist2.GetTopVolume() ) |
---|
119 | ) |
---|
120 | || |
---|
121 | ( |
---|
122 | ( hist1.GetTopVolume() == rh.hist2.GetTopVolume() ) && |
---|
123 | ( hist2.GetTopVolume() == rh.hist1.GetTopVolume() ) |
---|
124 | ) |
---|
125 | ) |
---|
126 | result = true; |
---|
127 | return result; |
---|
128 | } |
---|
129 | |
---|
130 | |
---|
131 | OlapEventAction::OlapEventAction(OlapRunAction * aRunAction) |
---|
132 | : theRunAction(aRunAction), dontDelete(false) |
---|
133 | { |
---|
134 | const G4VUserDetectorConstruction * aDet = |
---|
135 | G4RunManager::GetRunManager()->GetUserDetectorConstruction(); |
---|
136 | const OlapDetConstr * aConstDet = |
---|
137 | dynamic_cast<const OlapDetConstr *>(aDet); |
---|
138 | if (!aDet) { |
---|
139 | G4cerr << "OlapEventAction just can be used together with an OlapDetConstr!" |
---|
140 | << G4endl << "exiting ..." << G4endl; |
---|
141 | G4Exception("ERROR - OlapEventAction::OlapEventAction()"); |
---|
142 | } |
---|
143 | theDet = const_cast<OlapDetConstr *>(aConstDet); |
---|
144 | } |
---|
145 | |
---|
146 | |
---|
147 | OlapEventAction::~OlapEventAction() |
---|
148 | { |
---|
149 | } |
---|
150 | |
---|
151 | |
---|
152 | void OlapEventAction::BeginOfEventAction(const G4Event*) |
---|
153 | { |
---|
154 | while (!ABSteps.empty()) |
---|
155 | { |
---|
156 | if (! dontDelete ) |
---|
157 | delete ABSteps.back(); |
---|
158 | ABSteps.pop_back(); |
---|
159 | } |
---|
160 | |
---|
161 | while (!BASteps.empty()) |
---|
162 | { |
---|
163 | if (! dontDelete ) |
---|
164 | delete BASteps.back(); |
---|
165 | BASteps.pop_back(); |
---|
166 | } |
---|
167 | |
---|
168 | dontDelete=false; |
---|
169 | } |
---|
170 | |
---|
171 | |
---|
172 | void OlapEventAction::EndOfEventAction(const G4Event* anEvent) |
---|
173 | { |
---|
174 | |
---|
175 | const OlapGenerator * aGen = dynamic_cast<const OlapGenerator *> |
---|
176 | (G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction()); |
---|
177 | OlapGenerator * aGenerator=0; |
---|
178 | if ( aGen ) |
---|
179 | { |
---|
180 | aGenerator = const_cast<OlapGenerator*>(aGen); |
---|
181 | } |
---|
182 | else |
---|
183 | { |
---|
184 | G4cerr << "Warning: Primary generator is not OlapGenerator!" << G4endl |
---|
185 | << " Overlap Detection will not work!" << G4endl; |
---|
186 | return; |
---|
187 | } |
---|
188 | |
---|
189 | if (! (anEvent->GetEventID() % 1000)) |
---|
190 | { |
---|
191 | G4cout << "Event=" << anEvent->GetEventID() << " : " |
---|
192 | << theDet->GetNewWorld()->GetLogicalVolume()->GetDaughter(0)->GetName() |
---|
193 | << G4endl; |
---|
194 | } |
---|
195 | |
---|
196 | // try the starting points of AB and BA |
---|
197 | // they must be in the theWorldLV volume ... |
---|
198 | // G4LogicalVolume * startLV = |
---|
199 | // (ABSteps[0])->theHist.GetTopVolume()->GetLogicalVolume(); |
---|
200 | |
---|
201 | G4int axx = -1; |
---|
202 | if (aGenerator) |
---|
203 | axx = aGenerator->GetAxis(); |
---|
204 | |
---|
205 | G4double distAB, distBA; |
---|
206 | distAB = (ABSteps[0])->thePoint[axx] |
---|
207 | - (ABSteps[ABSteps.size()-1])->thePoint[axx]; |
---|
208 | distBA = (BASteps[0])->thePoint[axx] |
---|
209 | - (BASteps[BASteps.size()-1])->thePoint[axx]; |
---|
210 | #ifdef OLAP_DEBUG |
---|
211 | G4cout << "OlapEventAction::EndOfEventAction(): " |
---|
212 | << "deltaAB=" |
---|
213 | << (ABSteps[0])->thePoint[axx] |
---|
214 | - (ABSteps[ABSteps.size()-1])->thePoint[axx] |
---|
215 | << " deltaBA=" |
---|
216 | << (BASteps[0])->thePoint[axx] |
---|
217 | - (BASteps[BASteps.size()-1])->thePoint[axx] |
---|
218 | << G4endl; |
---|
219 | #endif |
---|
220 | |
---|
221 | G4double tolerance = OlapManager::GetOlapManager()->Delta(); |
---|
222 | |
---|
223 | if ( (ABSteps[ABSteps.size()-1])->theHist.GetTopVolume()->GetLogicalVolume() != |
---|
224 | theDet->GetNewWorld()->GetLogicalVolume() || |
---|
225 | (BASteps[BASteps.size()-1])->theHist.GetTopVolume()->GetLogicalVolume() != |
---|
226 | theDet->GetNewWorld()->GetLogicalVolume() |
---|
227 | ) |
---|
228 | { |
---|
229 | #ifdef OLAP_DEBUG_2 |
---|
230 | G4cerr << "=============overlap: daughter protrudes mother" << G4endl; |
---|
231 | #endif |
---|
232 | //G4cerr << ABSteps[0]->theHist << G4endl; |
---|
233 | OlapInfo * oi = |
---|
234 | new OlapInfo( ABSteps[ABSteps.size()-1]->theHist, |
---|
235 | BASteps[BASteps.size()-1]->theHist, |
---|
236 | ABSteps[ABSteps.size()-1]->thePoint, |
---|
237 | BASteps[BASteps.size()-1]->thePoint, |
---|
238 | axx, |
---|
239 | OlapManager::GetOlapManager()->GetOriginalWorld() ); |
---|
240 | oi->info = "daughter protrudes mother"; |
---|
241 | oi->probNot = false; // it's an overlap |
---|
242 | oi->stAB = ABSteps; |
---|
243 | oi->stBA = BASteps; |
---|
244 | dontDelete = true; |
---|
245 | |
---|
246 | if ( !InsertOI(oi) ) |
---|
247 | { // already detected! |
---|
248 | delete oi; |
---|
249 | } |
---|
250 | return; |
---|
251 | } |
---|
252 | |
---|
253 | // just in case, not the OlapGenerator is used but |
---|
254 | // another particle generator |
---|
255 | if (!BASteps.size()) |
---|
256 | return; |
---|
257 | |
---|
258 | // do we have as many steps from AB as from BA? |
---|
259 | if (ABSteps.size() != BASteps.size()) { |
---|
260 | #ifdef OLAP_DEBUG |
---|
261 | G4cerr << "overlap: different no of steps AB, BA" << G4endl; |
---|
262 | #endif |
---|
263 | //G4cerr << ABSteps[0]->theHist << G4endl; |
---|
264 | //G4cerr << "in that case we're are not logging (yet) ..." << G4endl; |
---|
265 | OlapInfo * oi = |
---|
266 | new OlapInfo( ABSteps[1]->theHist, |
---|
267 | BASteps[1]->theHist, |
---|
268 | ABSteps[1]->thePoint, |
---|
269 | BASteps[1]->thePoint, |
---|
270 | axx, |
---|
271 | OlapManager::GetOlapManager()->GetOriginalWorld() ); |
---|
272 | // copy all steps and prohibit deletion of their pointers ... |
---|
273 | oi->stAB = ABSteps; |
---|
274 | oi->stBA = BASteps; |
---|
275 | oi->probNot = true; // probably not an overlap, just numerics .... |
---|
276 | dontDelete = true; |
---|
277 | |
---|
278 | std::ostringstream os; |
---|
279 | os << "A,B diff. steps AB:" << ABSteps.size() |
---|
280 | << " BA:" << BASteps.size() << '\0'; |
---|
281 | |
---|
282 | oi->info = G4String(os.str()); |
---|
283 | if ( !InsertOI(oi) ) |
---|
284 | { |
---|
285 | delete oi; |
---|
286 | } |
---|
287 | return; |
---|
288 | } |
---|
289 | |
---|
290 | // now forget about the first and the last points of AB and BA & |
---|
291 | // compare the remaining points of AB with the remaining points |
---|
292 | // of BA where the points of BA are in descending order while |
---|
293 | // the points of BA are in ascending order |
---|
294 | G4int siz = ABSteps.size(); |
---|
295 | G4int j = siz-3; |
---|
296 | for (G4int i = 1; i<(siz-1); i++) |
---|
297 | { |
---|
298 | G4double delta = |
---|
299 | (ABSteps[i]->thePoint - BASteps[j]->thePoint).mag() ; |
---|
300 | if (delta>tolerance) |
---|
301 | { |
---|
302 | OlapInfo * oi = |
---|
303 | new OlapInfo( ABSteps[i]->theHist, |
---|
304 | BASteps[j]->theHist, |
---|
305 | ABSteps[i]->thePoint, |
---|
306 | BASteps[j]->thePoint, |
---|
307 | axx, |
---|
308 | OlapManager::GetOlapManager()->GetOriginalWorld() ); |
---|
309 | #ifdef OLAP_DEBUG |
---|
310 | G4cerr << "overlap: delta=" << delta << G4endl; |
---|
311 | G4cerr << *oi << G4endl; |
---|
312 | #endif |
---|
313 | |
---|
314 | if ( !InsertOI(oi) ) |
---|
315 | { // already detected! |
---|
316 | delete oi; |
---|
317 | } |
---|
318 | |
---|
319 | /* |
---|
320 | G4cerr << "overlap: delta=" << delta << G4endl; |
---|
321 | G4cerr << *(ABSteps[i]) |
---|
322 | << *(BASteps[j]); |
---|
323 | */ |
---|
324 | } |
---|
325 | j--; |
---|
326 | } |
---|
327 | |
---|
328 | // G4cout << "From A to B:" << G4endl; |
---|
329 | // G4cout << ABSteps << G4endl; |
---|
330 | // G4cout << "From B to A:" << G4endl; |
---|
331 | // G4cout << BASteps << G4endl; |
---|
332 | } |
---|
333 | |
---|
334 | G4bool OlapEventAction::InsertOI(OlapInfo * oi) |
---|
335 | { |
---|
336 | std::vector<OlapInfo*>::iterator it = theRunAction->theOlaps.begin(); |
---|
337 | G4bool flag(false); |
---|
338 | while (it != theRunAction->theOlaps.end()) |
---|
339 | { |
---|
340 | if (*oi== *(*it)) { // already detected overlap |
---|
341 | flag=true; |
---|
342 | if ((*it)->probNot != oi->probNot ) |
---|
343 | { // probably a real overlap, not numerics |
---|
344 | flag=false; |
---|
345 | (*it)->probNot = false; |
---|
346 | oi->probNot = false; |
---|
347 | } |
---|
348 | break; |
---|
349 | } |
---|
350 | it++; |
---|
351 | } |
---|
352 | |
---|
353 | if (flag) |
---|
354 | { |
---|
355 | #ifdef OLAP_DEBUG |
---|
356 | G4cout << "===================================" << G4endl; |
---|
357 | G4cout << "No of detected Overlaps: " << theRunAction->theOlaps.size() |
---|
358 | << G4endl; |
---|
359 | G4cout << "===================================" << G4endl; |
---|
360 | #endif |
---|
361 | return false; // already detected overlap |
---|
362 | } |
---|
363 | |
---|
364 | theRunAction->theOlaps.push_back(oi); |
---|
365 | |
---|
366 | #ifdef OLAP_DEBUG |
---|
367 | G4cout << "===================================" << G4endl; |
---|
368 | G4cout << "No of detected Overlaps: " << theRunAction->theOlaps.size() |
---|
369 | << G4endl; |
---|
370 | G4cout << "===================================" << G4endl; |
---|
371 | G4cout << oi->hist1 << G4endl; |
---|
372 | G4cout << "CpNo: " << oi->hist1.GetTopVolume()->GetCopyNo() << G4endl; |
---|
373 | G4cout << "Pos: " << oi->hist1.GetTopTransform().NetTranslation() |
---|
374 | << G4endl; |
---|
375 | G4cout << "Rot: " << oi->hist1.GetTopTransform().NetRotation().phiX() |
---|
376 | << " " << oi->hist1.GetTopTransform().NetRotation().phiY() << " " |
---|
377 | << oi->hist1.GetTopTransform().NetRotation().phiZ() << " " |
---|
378 | << oi->hist1.GetTopTransform().NetRotation().thetaX() << " " |
---|
379 | << oi->hist1.GetTopTransform().NetRotation().thetaY() << " " |
---|
380 | << oi->hist1.GetTopTransform().NetRotation().thetaZ() |
---|
381 | << G4endl ; |
---|
382 | G4cout << oi->v1 << G4endl << "----" << G4endl; |
---|
383 | |
---|
384 | G4cout << oi->hist2 << G4endl; |
---|
385 | G4cout << "CpNo: " << oi->hist2.GetTopVolume()->GetCopyNo() << G4endl; |
---|
386 | G4cout << "Pos: " << oi->hist2.GetTopTransform().NetTranslation() |
---|
387 | << G4endl; |
---|
388 | G4cout << "Rot: " << oi->hist2.GetTopTransform().NetRotation().phiX() |
---|
389 | << " " << oi->hist2.GetTopTransform().NetRotation().phiY() << " " |
---|
390 | << oi->hist2.GetTopTransform().NetRotation().phiZ() << " " |
---|
391 | << oi->hist2.GetTopTransform().NetRotation().thetaX() << " " |
---|
392 | << oi->hist2.GetTopTransform().NetRotation().thetaY() << " " |
---|
393 | << oi->hist2.GetTopTransform().NetRotation().thetaZ() |
---|
394 | << G4endl ; |
---|
395 | G4cout << oi->v2 << G4endl << "----" << G4endl; |
---|
396 | |
---|
397 | #endif |
---|
398 | |
---|
399 | static char * c = getenv("OLAP_LOG_EVENT"); |
---|
400 | if (c) |
---|
401 | G4cerr << *oi << G4endl; |
---|
402 | |
---|
403 | return true; // new overlap detected! |
---|
404 | } |
---|