BOSS 7.0.6
BESIII Offline Software System
Loading...
Searching...
No Matches
G4SubtractionSolid.cc
Go to the documentation of this file.
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: G4SubtractionSolid.cc,v 1.2 2015/03/17 05:50:57 sunss Exp $
28// GEANT4 tag $Name: BesSim-00-01-24 $
29//
30// Implementation of methods for the class G4IntersectionSolid
31//
32// History:
33//
34// 14.10.98 V.Grichine: implementation of the first version
35// 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v)
36// 02.08.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
37// while -> do-while & surfaceA limitations
38// 13.09.00 V.Grichine: bug fixed in SurfaceNormal(p), p can be inside
39//
40// --------------------------------------------------------------------
41
42#include "G4SubtractionSolid.hh"
43
44#include "G4VoxelLimits.hh"
45#include "G4VPVParameterisation.hh"
46#include "G4GeometryTolerance.hh"
47
48#include "G4VGraphicsScene.hh"
49#include "G4Polyhedron.hh"
50#include "G4NURBS.hh"
51// #include "G4NURBSbox.hh"
52
53#include <sstream>
54
55///////////////////////////////////////////////////////////////////
56//
57// Transfer all data members to G4BooleanSolid which is responsible
58// for them. pName will be in turn sent to G4VSolid
59
60G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
61 G4VSolid* pSolidA ,
62 G4VSolid* pSolidB )
63 : G4BooleanSolid(pName,pSolidA,pSolidB)
64{
65}
66
67///////////////////////////////////////////////////////////////
68//
69// Constructor
70
71G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
72 G4VSolid* pSolidA ,
73 G4VSolid* pSolidB ,
74 G4RotationMatrix* rotMatrix,
75 const G4ThreeVector& transVector )
76 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
77{
78}
79
80///////////////////////////////////////////////////////////////
81//
82// Constructor
83
84G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
85 G4VSolid* pSolidA ,
86 G4VSolid* pSolidB ,
87 const G4Transform3D& transform )
88 : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
89{
90}
91
92//////////////////////////////////////////////////////////////////
93//
94// Fake default constructor - sets only member data and allocates memory
95// for usage restricted to object persistency.
96
97G4SubtractionSolid::G4SubtractionSolid( __void__& a )
98 : G4BooleanSolid(a)
99{
100}
101
102///////////////////////////////////////////////////////////////
103//
104// Destructor
105
106G4SubtractionSolid::~G4SubtractionSolid()
107{
108}
109
110///////////////////////////////////////////////////////////////
111//
112// CalculateExtent
113
114G4bool
115G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
116 const G4VoxelLimits& pVoxelLimit,
117 const G4AffineTransform& pTransform,
118 G4double& pMin,
119 G4double& pMax ) const
120{
121 // Since we cannot be sure how much the second solid subtracts
122 // from the first, we must use the first solid's extent!
123
124 return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
125 pTransform, pMin, pMax );
126}
127
128/////////////////////////////////////////////////////
129//
130// Touching ? Empty subtraction ?
131
132EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const
133{
134 EInside positionA = fPtrSolidA->Inside(p);
135 if (positionA == kOutside) return kOutside;
136
137 EInside positionB = fPtrSolidB->Inside(p);
138
139 if(positionA == kInside && positionB == kOutside)
140 {
141 return kInside ;
142 }
143 else
144 {
145 if(( positionA == kInside && positionB == kSurface) ||
146 ( positionB == kOutside && positionA == kSurface) ||
147 ( positionA == kSurface && positionB == kSurface &&
148 ( fPtrSolidA->SurfaceNormal(p) -
149 fPtrSolidB->SurfaceNormal(p) ).mag2() >
150 1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
151 {
152 return kSurface;
153 }
154 else
155 {
156 return kOutside;
157 }
158 }
159}
160
161//////////////////////////////////////////////////////////////
162//
163// SurfaceNormal
164
165G4ThreeVector
166G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const
167{
168 G4ThreeVector normal;
169 if( Inside(p) == kOutside )
170 {
171#ifdef G4BOOLDEBUG
172 G4cout << "WARNING - Invalid call [1] in "
173 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
174 << " Point p is inside !" << G4endl;
175 G4cout << " p = " << p << G4endl;
176 G4cerr << "WARNING - Invalid call [1] in "
177 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
178 << " Point p is inside !" << G4endl;
179 G4cerr << " p = " << p << G4endl;
180#endif
181 }
182 else
183 {
184 if( fPtrSolidA->Inside(p) == kSurface &&
185 fPtrSolidB->Inside(p) != kInside )
186 {
187 normal = fPtrSolidA->SurfaceNormal(p) ;
188 }
189 else if( fPtrSolidA->Inside(p) == kInside &&
190 fPtrSolidB->Inside(p) != kOutside )
191 {
192 normal = -fPtrSolidB->SurfaceNormal(p) ;
193 }
194 else
195 {
196 if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
197 {
198 normal = fPtrSolidA->SurfaceNormal(p) ;
199 }
200 else
201 {
202 normal = -fPtrSolidB->SurfaceNormal(p) ;
203 }
204#ifdef G4BOOLDEBUG
205 if(Inside(p) == kInside)
206 {
207 G4cout << "WARNING - Invalid call [2] in "
208 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
209 << " Point p is inside !" << G4endl;
210 G4cout << " p = " << p << G4endl;
211 G4cerr << "WARNING - Invalid call [2] in "
212 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
213 << " Point p is inside !" << G4endl;
214 G4cerr << " p = " << p << G4endl;
215 }
216#endif
217 }
218 }
219 return normal;
220}
221
222/////////////////////////////////////////////////////////////
223//
224// The same algorithm as in DistanceToIn(p)
225
226G4double
227G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p,
228 const G4ThreeVector& v ) const
229{
230 G4double dist = 0.0,disTmp = 0.0 ;
231
232#ifdef G4BOOLDEBUG
233 if( Inside(p) == kInside )
234 {
235 G4cout << "WARNING - Invalid call in "
236 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
237 << " Point p is inside !" << G4endl;
238 G4cout << " p = " << p << G4endl;
239 G4cout << " v = " << v << G4endl;
240 G4cerr << "WARNING - Invalid call in "
241 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
242 << " Point p is inside !" << G4endl;
243 G4cerr << " p = " << p << G4endl;
244 G4cerr << " v = " << v << G4endl;
245 }
246#endif
247
248 // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B
249 if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B
250 {
251 dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
252
253 if( fPtrSolidA->Inside(p+dist*v) != kInside )
254 {
255 do
256 {
257 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
258
259 if(disTmp == kInfinity)
260 {
261 return kInfinity ;
262 }
263 dist += disTmp ;
264
265 if( Inside(p+dist*v) == kOutside )
266 {
267 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
268 dist += disTmp ;
269 }
270 }
271 while( Inside(p+dist*v) == kOutside ) ;
272 }
273 }
274 else // p outside A, start in A
275 {
276 dist = fPtrSolidA->DistanceToIn(p,v) ;
277
278 if( dist == kInfinity ) // past A, hence past A\B
279 {
280 return kInfinity ;
281 }
282 else
283 {
284 int loopflag=0;
285 while( Inside(p+dist*v) == kOutside ) // pushing loop
286 {
287 loopflag++;
288 if(loopflag>10000) {
289 std::cout<<"G4SubtractionSolid::DistanceToIn loop >10000"<<std::endl;
290 return kInfinity ;
291 }
292 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
293 dist += disTmp ;
294
295 if( Inside(p+dist*v) == kOutside )
296 {
297 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
298
299 if(disTmp == kInfinity) // past A, hence past A\B
300 {
301 return kInfinity ;
302 }
303 dist += disTmp ;
304 }
305 }
306 }
307 }
308
309 return dist ;
310}
311
312////////////////////////////////////////////////////////
313//
314// Approximate nearest distance from the point p to the intersection of
315// two solids. It is usually underestimated from the point of view of
316// isotropic safety
317
318G4double
319G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const
320{
321 G4double dist=0.0;
322
323#ifdef G4BOOLDEBUG
324 if( Inside(p) == kInside )
325 {
326 G4cout << "WARNING - Invalid call in "
327 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
328 << " Point p is inside !" << G4endl;
329 G4cout << " p = " << p << G4endl;
330 G4cerr << "WARNING - Invalid call in "
331 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
332 << " Point p is inside !" << G4endl;
333 G4cerr << " p = " << p << G4endl;
334 }
335#endif
336
337 if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1
338 ( fPtrSolidB->Inside(p) != kOutside) )
339 {
340 dist= fPtrSolidB->DistanceToOut(p) ;
341 }
342 else
343 {
344 dist= fPtrSolidA->DistanceToIn(p) ;
345 }
346
347 return dist;
348}
349
350//////////////////////////////////////////////////////////
351//
352// The same algorithm as DistanceToOut(p)
353
354G4double
355G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p,
356 const G4ThreeVector& v,
357 const G4bool calcNorm,
358 G4bool *validNorm,
359 G4ThreeVector *n ) const
360{
361#ifdef G4BOOLDEBUG
362 if( Inside(p) == kOutside )
363 {
364 G4cout << "Position:" << G4endl << G4endl;
365 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
366 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
367 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
368 G4cout << "Direction:" << G4endl << G4endl;
369 G4cout << "v.x() = " << v.x() << G4endl;
370 G4cout << "v.y() = " << v.y() << G4endl;
371 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
372 G4cout << "WARNING - Invalid call in "
373 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
374 << " Point p is outside !" << G4endl;
375 G4cout << " p = " << p << G4endl;
376 G4cout << " v = " << v << G4endl;
377 G4cerr << "WARNING - Invalid call in "
378 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
379 << " Point p is outside !" << G4endl;
380 G4cerr << " p = " << p << G4endl;
381 G4cerr << " v = " << v << G4endl;
382 }
383#endif
384
385 G4double distout;
386 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
387 G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
388 if(distB < distA)
389 {
390 if(calcNorm)
391 {
392 *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
393 *validNorm = false ;
394 }
395 distout= distB ;
396 }
397 else
398 {
399 distout= distA ;
400 }
401 return distout;
402}
403
404//////////////////////////////////////////////////////////////
405//
406// Inverted algorithm of DistanceToIn(p)
407
408G4double
409G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const
410{
411 G4double dist=0.0;
412
413 if( Inside(p) == kOutside )
414 {
415#ifdef G4BOOLDEBUG
416 G4cout << "WARNING - Invalid call in "
417 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
418 << " Point p is outside" << G4endl;
419 G4cout << " p = " << p << G4endl;
420 G4cerr << "WARNING - Invalid call in "
421 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
422 << " Point p is outside" << G4endl;
423 G4cerr << " p = " << p << G4endl;
424#endif
425 }
426 else
427 {
428 dist= std::min(fPtrSolidA->DistanceToOut(p),
429 fPtrSolidB->DistanceToIn(p) ) ;
430 }
431 return dist;
432}
433
434//////////////////////////////////////////////////////////////
435//
436//
437
438G4GeometryType G4SubtractionSolid::GetEntityType() const
439{
440 return G4String("G4SubtractionSolid");
441}
442
443//////////////////////////////////////////////////////////////
444//
445//
446
447void
448G4SubtractionSolid::ComputeDimensions( G4VPVParameterisation*,
449 const G4int,
450 const G4VPhysicalVolume* )
451{
452}
453
454/////////////////////////////////////////////////
455//
456//
457
458void
459G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
460{
461 scene.AddSolid (*this);
462}
463
464////////////////////////////////////////////////////
465//
466//
467
468G4Polyhedron*
469G4SubtractionSolid::CreatePolyhedron () const
470{
471 G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
472 G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
473 if (pA && pB)
474 {
475 G4Polyhedron* resultant = new G4Polyhedron (pA->subtract(*pB));
476 return resultant;
477 }
478 else
479 {
480 std::ostringstream oss;
481 oss << "Solid - " << GetName()
482 << " - one of the Boolean components has no" << G4endl
483 << " corresponding polyhedron. Returning NULL !";
484 G4Exception("G4SubtractionSolid::CreatePolyhedron()", "InvalidSetup",
485 JustWarning, oss.str().c_str());
486 return 0;
487 }
488}
489
490/////////////////////////////////////////////////////////
491//
492//
493
494G4NURBS*
495G4SubtractionSolid::CreateNURBS () const
496{
497 // Take into account boolean operation - see CreatePolyhedron.
498 // return new G4NURBSbox (1.0, 1.0, 1.0);
499 return 0;
500}
pMin
const double pMin
Definition: CheckRung/check/dEdxVsP.cxx:24
pMax
const double pMax
Definition: CheckRung/check/dEdxVsP.cxx:25
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition: KarLud.h:35