source: trunk/examples/extended/geometry/olap/src/OlapEventAction.cc@ 1010

Last change on this file since 1010 was 807, checked in by garnier, 17 years ago

update
File size: 13.5 KB
RevLine 
[807]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
59std::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
83OlapInfo::~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
99G4bool 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
131OlapEventAction::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
147OlapEventAction::~OlapEventAction()
148{
149}
150
151
152void 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
172void 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
334G4bool 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}
Note: See TracBrowser for help on using the repository browser.