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

Last change on this file since 1279 was 1230, checked in by garnier, 14 years ago

update to geant4.9.3
File size: 13.6 KB
Line 
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: geant4-09-03-cand-01 $
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.