source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/G4Detector/G4fresnellens/src/ShapeOfLens.cc @ 117

Last change on this file since 117 was 117, checked in by moretto, 11 years ago

ESAF version compilable on mac OS
File size: 19.0 KB
Line 
1#include "ShapeOfLens.h"
2#include "FresnelSurface.h"
3#include "LensUtils.h"
4#include "FresnelOutput_def.h"
5
6#include <globals.hh>
7#include <G4VoxelLimits.hh>
8#include <G4AffineTransform.hh>
9#include <G4GeometryTolerance.hh>
10#include <G4OpBoundaryProcess.hh>
11#include <G4Tubs.hh>
12#include <G4EventManager.hh>
13#include <G4Track.hh>
14#include <Randomize.hh>
15
16using namespace G4FresnelLens;
17using namespace CLHEP;
18using std::min;
19using std::max;
20using std::swap;
21
22//#define DEBUG_WNF
23//#define DEBUG
24//#define DEBUG_SLANTWALLS
25
26///////////////////////////////////////////////////////////////////////////////
27//
28// constructor - check parameters
29//__________________________________________________________________________________
30ShapeOfLens::ShapeOfLens(const G4String& pName,
31                               G4double pDz,
32                               G4double pR,
33                               SSurface* upsurf, SSurface* downsurf)
34: G4VSolid(pName),
35__INHERITCONSTRUCTOR("ShapeOfLens")
36fCubicVolume(0.), //fpPolyhedron(0),
37fUpSurface(upsurf), fDownSurface(downsurf)
38#ifdef SLANTWALLS
39,isSlant(false), theSlantSigma(0.)
40#endif
41{
42    top_volume = new G4Tubs("top_volume",0.,pR,pDz,0.,360.);
43    kTolerance=G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
44    kHalfTolerance=kTolerance*.5;
45    fR=min(fUpSurface->GetR(), fDownSurface->GetR());
46    fZ1=fDownSurface->GetMinZ();
47    fZ2=fUpSurface->GetMaxZ();
48    fCylinder.SetR(fR);
49    fCylinder.SetZ(fDownSurface->GetEdgeZ(), fUpSurface->GetEdgeZ());
50
51    theLastPosition.set(DBL_MAX, DBL_MAX, DBL_MAX);
52    theLastInside=kOutside;
53
54    fSurfaceArea = 0.;
55
56    __REG_ACTION("ShapeOfLens");
57}
58
59//__________________________________________________________________________________
60ShapeOfLens::ShapeOfLens( __void__& a ) : G4VSolid(a),  fCubicVolume(0) {
61    fSurfaceArea = 0.;
62}
63
64//__________________________________________________________________________________
65ShapeOfLens::~ShapeOfLens() {
66    // top_volume is deleted by Geant4
67    //delete top_volume;
68    if ( fUpSurface ) delete fUpSurface;
69    if ( fDownSurface ) delete fDownSurface;
70}
71
72//__________________________________________________________________________________
73EInside ShapeOfLens::Inside(const G4ThreeVector& p, int* surface, FSIntersectedSegment* intSeg) const {
74    #ifdef DEBUG
75    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
76    printf("Inside called (%s)\n", GetName().data());
77    printf("din pos(%f, %f, %f)\n",p.x(),p.y(),p.z());
78    printf("rho=%f z=%f\n",p.perp(),p.z());
79    #endif
80    //if ( VectorIsZero(p) ) printf("zero\n");
81
82    if ( CompareVectors(p, theLastPosition, kHalfTolerance) ) {
83        if ( intSeg ) *intSeg=theLastInsideSegment;
84        #ifdef DEBUG
85        printf("vector is cached.\n");
86        printf("return %s\n", theLastInside==kInside?"inside":(theLastInside==kOutside?"out":"surf"));
87        #endif
88        return theLastInside;
89    }
90    theLastInsideSegment.Reset();
91    theLastPosition=p;
92
93    double z = p.z();
94    if ( z<fZ1-kHalfTolerance || z>fZ2+kHalfTolerance ) {
95        #ifdef DEBUG
96        printf("return out, for sure\n");
97        #endif
98        return theLastInside=kOutside;
99    }
100
101    double rho = p.perp();
102    if ( rho > fR+kHalfTolerance ) {
103        #ifdef DEBUG
104        printf("return out, on side\n");
105        #endif
106        return theLastInside=kOutside;
107    }
108
109
110    EInside ins=fCylinder.Inside(rho, z, &theLastInsideSegment);
111    if (ins==kSurface) {
112        if (surface) *surface = 0;
113        #ifdef DEBUG
114        printf("return surf, cyl\n");
115        #endif
116        return theLastInside=kSurface;
117    }
118
119    EInside in = fDownSurface->Inside(rho, z, &theLastInsideSegment);
120    if ( in==kInside ) {
121        #ifdef DEBUG
122        printf("return out, down\n");
123        #endif
124        return theLastInside=kOutside;
125    }
126    if ( in==kSurface) {
127        if ( intSeg ) {
128            *intSeg=theLastInsideSegment;
129            if (intSeg->surf) printf("cyl err\n");
130        }
131        if (surface) *surface = 1;
132        #ifdef DEBUG
133        printf("return surface, down\n");
134        #endif
135        return theLastInside=kSurface;
136    }
137
138    in = fUpSurface->Inside(rho, z, &theLastInsideSegment);
139    if ( intSeg && in==kSurface ) *intSeg=theLastInsideSegment;
140    if (surface) *surface = 2;
141    #ifdef DEBUG
142    printf("return %s, up\n", in==kInside?"inside":(in==kOutside?"out":"in"));
143    #endif
144    return theLastInside=in;
145}
146
147//__________________________________________________________________________________
148G4double ShapeOfLens::DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& v) const {
149    #ifdef DEBUG
150    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
151    printf("DistanceToIn(vec) called (%s)\n", GetName().data());
152    printf("din pos(%.10g, %.10g, %.10g) dir(%.10g, %.10g, %.10g)\n",p.x(),p.y(),p.z(),v.x(),v.y(),v.z());
153    printf("rho=%g z=%g, theta=%g phi=%g\n",p.perp(),p.z(), v.theta()/deg, v.phi()/deg);
154    #endif
155
156    double zp=p.z();
157    double vz=v.z();
158
159    theLastPositionOnSurface.set(kInfinity, kInfinity, kInfinity);
160    theLastSurface=0;
161    if ( (zp>=fZ2 && vz>=0.) || (zp<=fZ1 && vz<=0.) ) return kInfinity;
162
163    double perp=p.perp();
164    double smallSubst=0;
165    FSIntersectedSegment intUp, intDown;
166    theLastPositionOnSurface.set(kInfinity, kInfinity, kInfinity);
167    if ( perp<(fR-kHalfTolerance) ) {
168        if ( zp>=(fUpSurface->GetBottomZ(perp)-kHalfTolerance) ) {
169            EInside in = fUpSurface->Inside(perp, zp);
170            #ifdef DEBUG
171            printf("in_up: %s\n", in==kSurface?"sur":( in==kOutside?"out":"in"));
172            #endif
173            if ( in!=kInside) {
174                G4ThreeVector p1=p;
175                if ( in==kSurface ) p1-=v*kHalfTolerance;
176                double dd=fUpSurface->FindNearestIntersection(p1, v, &intUp);
177                if ( in==kSurface ) dd-=kHalfTolerance;
178                if( dd<=kHalfTolerance ) dd=smallSubst;;
179                theLastPositionOnSurface=p+v*dd;
180                theLastSurfaceSegment=intUp;
181                theLastSurface=2;
182                #ifdef DEBUG
183                printf("dist up=%f\n",dd);
184                if ( dd<kInfinity ) {
185                    printf("next point (%f, %f, %f)\n",theLastPositionOnSurface.x(), theLastPositionOnSurface.y(), theLastPositionOnSurface.z());
186                }
187                #endif
188                return dd;
189            }
190        }
191        if ( zp<=(fDownSurface->GetUpperZ(perp)+kHalfTolerance) ) {
192            EInside in = fDownSurface->Inside(perp, zp);
193            #ifdef DEBUG
194            printf("in_down: %s\n", in==kSurface?"sur":( in==kOutside?"out":"in"));
195            #endif
196            if ( in!=kOutside) {
197                G4ThreeVector p1=p;
198                if ( in==kSurface ) p1-=v*kHalfTolerance;
199                double dd=fDownSurface->FindNearestIntersection(p, v, &intDown);
200                if ( in==kSurface ) dd-=kHalfTolerance;
201                if( dd<=kHalfTolerance ) dd=smallSubst;
202                theLastPositionOnSurface=p+v*dd;
203                theLastSurfaceSegment=intDown;
204                theLastSurface=1;
205                #ifdef DEBUG
206                printf("dist down=%f\n",dd);
207                if ( dd<kInfinity ) {
208                    printf("next point (%f, %f, %f)\n",theLastPositionOnSurface.x(), theLastPositionOnSurface.y(), theLastPositionOnSurface.z());
209                }
210                #endif
211                return dd;
212            }
213        }
214
215        Warning_("Point is inside, but calling DistanceToIn.");
216        int id=G4EventManager::GetEventManager()->
217                               GetTrackingManager()->
218                               GetTrack()->GetTrackID();
219        printf("Event %i, %s\n", id, GetName().data());
220        printf("pos(%.10g,%.10g,%.10g)\n", p.x(), p.y(), p.z());
221        printf("dif(%.10g,%.10g,%.10g)\n", v.x(), v.y(), v.z());
222
223        return kInfinity;
224    }
225
226    double cosalpha = p.x()*v.x()+p.y()*v.y();
227    if ( cosalpha>=.0 ) return kInfinity;
228    // calculate distance to cylinder
229    FSIntersectedSegment intCyl;
230    double dist_cyl  = fCylinder.FindNearestIntersection(p, v, &intCyl);
231    double dist_up   = fUpSurface->FindNearestIntersection(p, v, &intUp);
232    double dist_down = fDownSurface->FindNearestIntersection(p, v, &intDown);
233
234    #ifdef DEBUG
235    printf("dc: %e \tdu: %e \t dd:%e\n",dist_cyl,dist_up,dist_down);
236    #endif
237
238    double mindist;
239    if ( dist_up<dist_down ) {
240        if ( dist_up<dist_cyl ){
241            mindist=dist_up;
242            theLastSurfaceSegment=intUp;
243            theLastSurface=2;
244        }
245        else {
246            mindist=dist_cyl;
247            theLastSurfaceSegment=intCyl;
248            theLastSurface=0;
249        }
250    }
251    else {
252        if ( dist_down<dist_cyl ){
253            mindist=dist_down;
254            theLastSurfaceSegment=intDown;
255            theLastSurface=1;
256        }
257        else {
258            mindist=dist_cyl;
259            theLastSurfaceSegment=intCyl;
260            theLastSurface=0;
261        }
262    }
263    mindist=mindist>kTolerance?dist_cyl:smallSubst;
264    theLastPositionOnSurface=p+v*mindist;
265    return mindist;
266}
267
268//__________________________________________________________________________________
269G4double ShapeOfLens::DistanceToIn(const G4ThreeVector& p) const{
270    #ifdef DEBUG
271    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
272    printf("DistanceToIn(pos) called (%s)\n", GetName().data());
273    printf("din pos(%f, %f, %f)\n",p.x(),p.y(),p.z());
274    printf("rho=%f z=%f\n",p.perp(),p.z());
275    #endif
276    double zp=p.z();
277    if ( zp<=fCylinder.GetMinZ() ) return fDownSurface->DistanceToSurfSafeDown(p);
278    if ( zp>=fCylinder.GetMaxZ() ) return fUpSurface->DistanceToSurfSafeUp(p);
279    return fabs(p.perp()-fR);
280}
281
282//__________________________________________________________________________________
283G4double ShapeOfLens::DistanceToOut(const G4ThreeVector& p,
284                           const G4ThreeVector& v,
285                           const G4bool calcNorm,
286                                 G4bool *validNorm,
287                                 G4ThreeVector *n) const {
288    #ifdef DEBUG
289    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
290    printf("DistanceToOut(vec) called (%s)\n", GetName().data());
291    printf("dout (%f, %f, %f) (%f, %f, %f)\n",p.x(),p.y(),p.z(),v.x(),v.y(),v.z());    printf("rho=%f z=%f\n",p.perp(),p.z());
292    #endif
293    if ( calcNorm ) *validNorm = true;
294
295    G4ThreeVector p1=p;
296    double dist_shift(0.);
297    double perp=p.perp();
298    EInside in_up = fUpSurface->Inside(perp,p.z());
299    EInside in_down = fDownSurface->Inside(perp,p.z());
300    #ifdef DEBUG
301    printf("in_up: %s\n",in_up==kSurface?"sur":(in_up==kOutside?"out":"in"));
302    printf("in_down: %s\n",in_down==kSurface?"sur":(in_down==kOutside?"out":"in"));
303    #endif
304    if ( in_up==kSurface || in_down==kSurface ) {
305        #ifdef DEBUG
306        printf("point is on surface, shift it by kHalfTolerance\n");
307        #endif
308        p1+=kHalfTolerance*v;
309        dist_shift+=kHalfTolerance;
310    }
311
312    #ifdef DEBUG
313    printf("asking up\n");
314    #endif
315    FSIntersectedSegment intsegUp, intsegDown;
316    double dist_up = fUpSurface->FindNearestIntersection(p1, v, &intsegUp);
317    #ifdef DEBUG
318    printf("asking down\n");
319    #endif
320    double dist_down = fDownSurface->FindNearestIntersection(p1, v, &intsegDown);
321
322    #ifdef DEBUG
323      printf("du=%e dd=%e\n",dist_up,dist_down);
324      double du=dist_up+dist_shift;
325      double dd=dist_down+dist_shift;
326      printf("du(s)=%e dd=%e\n",du,dd);
327      if ( du<kInfinity ) {
328          G4ThreeVector newp=p+v*du;
329          printf("next point (up) (%f, %f, %f)\n",newp.x(), newp.y(), newp.z());
330      }
331      if ( dd<kInfinity ) {
332          G4ThreeVector newp=p+v*dd;
333          printf("next point (down) (%f, %f, %f)\n",newp.x(), newp.y(), newp.z());
334      }
335    #endif
336
337    if ( dist_up<kInfinity || dist_down<kInfinity ) {
338        if ( dist_up<dist_down ) {
339            if ( calcNorm ) {
340                intsegUp.Normal(p1+v*dist_up, +1., n);
341                #ifdef SLANTWALLS
342                if ( isSlant && n->z()==0. ) CorrectCylinderNormal(n);
343                #endif
344                #ifdef DEBUG
345                printf("normal is (%f, %f, %f)\n",n->x(), n->y(), n->z());
346                #endif
347            }
348            dist_shift+=dist_up;
349            return dist_shift<kHalfTolerance?kHalfTolerance:dist_shift+kHalfTolerance*0.5;
350        }
351
352        if ( dist_down<kInfinity ) {
353            if ( calcNorm ) {
354                intsegDown.Normal(p1+v*dist_down, -1., n);
355                #ifdef SLANTWALLS
356                if ( isSlant && n->z()==0.) CorrectCylinderNormal(n);
357                #endif
358                #ifdef DEBUG
359                printf("normal is (%f, %f, %f)\n",n->x(), n->y(), n->z());
360                #endif
361            }
362            dist_shift+=dist_down;
363            return dist_shift<kHalfTolerance?kHalfTolerance:dist_shift+kHalfTolerance*0.5;
364        }
365    }
366
367    double dist;
368    G4ThreeVector pp;
369    if ( calcNorm ) {
370        FSIntersectedSegment intCyl;
371        dist=fCylinder.FindNearestIntersection(p1, v, &intCyl);
372        pp=p1+v*dist;
373        intCyl.Normal(pp, -1., n);
374        #ifdef SLANTWALLS
375        if ( isSlant ) CorrectCylinderNormal(n);
376        #endif
377    }
378    else {
379        dist=fCylinder.FindNearestIntersection(p1, v);
380        pp=p1+v*dist;
381    }
382
383    if ( dist==kInfinity ) {
384        if ( dist_shift ) return dist_shift<kHalfTolerance?kHalfTolerance:dist_shift;
385        else {
386            #ifdef DEBUG_WNF
387                printf("shift is %g\n", dist_shift);
388                if ( fUpSurface->getDebug() ) exit(1);
389                fUpSurface->setDebug();
390                fDownSurface->setDebug();
391                DistanceToOut(p, v, calcNorm, validNorm, n);
392            #endif
393            Warning_("Way out not found.");
394            int id=G4EventManager::GetEventManager()->
395                                   GetTrackingManager()->
396                                   GetTrack()->GetTrackID();
397            printf("Event %i, %s\n", id, GetName().data());
398            printf("pos(%.10g,%.10g,%.10g)\n", p.x(), p.y(), p.z());
399            printf("dif(%.10g,%.10g,%.10g)\n", v.x(), v.y(), v.z());
400        }
401    }
402    dist+=dist_shift;
403    return dist<kHalfTolerance?kHalfTolerance:dist+kHalfTolerance*0.5;
404}
405
406//__________________________________________________________________________________
407G4double ShapeOfLens::DistanceToOut(const G4ThreeVector& p) const{
408    #ifdef DEBUG
409    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
410    printf("DistanceToOut(pos) called (%s)\n", GetName().data());
411    printf("din pos(%f, %f, %f)\n",p.x(),p.y(),p.z());
412    printf("rho=%f z=%f\n",p.perp(),p.z());
413    #endif
414    double d1=fDownSurface->DistanceToSurfSafeUp(p);
415    if ( d1<kTolerance ){
416        #ifdef DEBUG
417            printf("return dist to up: %g\n",d1);
418        #endif
419        return d1;
420    }
421    double d2=fUpSurface->DistanceToSurfSafeDown(p);
422    if ( d2<kTolerance ){
423        #ifdef DEBUG
424            printf("dist to up: %g\n", d1);
425            printf("return dist to down: %g\n",d2);
426        #endif
427        return d2;
428    }
429    double d3=fabs(fR-p.perp());
430
431    #ifdef DEBUG
432        printf("dist to up: %g\n", d1);
433        printf("dist to down: %g\n", d2);
434        printf("dist to cyl: %g\n", d3);
435        printf("return: %g\n", min(min(d1,d2),d3));
436    #endif
437    return min(min(d1,d2),d3);
438}
439
440//#define DEBUG_SAVENORMAL
441//__________________________________________________________________________________
442G4ThreeVector ShapeOfLens::SurfaceNormal( const G4ThreeVector& p ) const {
443    #ifdef DEBUG
444    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
445    printf("SurfaceNormal called (%s)\n", GetName().data());
446    printf("p (%f, %f)\n",p.perp(),p.z());
447    #endif
448    int surface(theLastSurface);
449    FSIntersectedSegment intseg;
450    G4ThreeVector norm(1.,0.,0.);
451    if ( CompareVectors(p, theLastPositionOnSurface, kHalfTolerance) ){
452        intseg=theLastSurfaceSegment;
453        #ifdef DEBUG
454        printf("Point is saved, using the last segment normal.\n");
455        #endif
456    }
457    else {
458        EInside ins=this->Inside(p, &surface, &intseg);
459        if ( ins!=kSurface ){
460            Warning_("Calculating SurfaceNormal not on surface.");
461            int id=G4EventManager::GetEventManager()->
462                                   GetTrackingManager()->
463                                   GetTrack()->GetTrackID();
464            printf("Event %i, %s\n", id, GetName().data());
465            printf("pos(%.10g,%.10g,%.10g)\n", p.x(), p.y(), p.z());
466            double rho=p.perp();
467            printf("dz1 %g, dz2 %g\n", p.z()-fUpSurface->PointOnSurf(rho),p.z()-fDownSurface->PointOnSurf(rho));
468            return norm;
469        }
470
471        #ifdef DEBUG
472        printf("Point is not saved. Dist is %g\n", (p-theLastPositionOnSurface).mag() );
473        #endif
474    }
475    #ifdef DEBUG
476    printf("    Point: %g %g %g\n",p.x(), p.y(), p.z());
477    printf("    Saved point: %g %g %g\n",theLastPositionOnSurface.x(), theLastPositionOnSurface.y(), theLastPositionOnSurface.z());
478    printf("    Saved segment: type %2i, r(%g, %g), z(%g, %g)\n", intseg.type, intseg.r1, intseg.r2, intseg.z1, intseg.z2);
479    #endif
480    intseg.Normal(p, surface<2?-1.:1., &norm);
481    #ifdef DEBUG
482    double  mag=norm.mag();
483    printf("Returned normal: (%g, %g, %g) (%g, %g) len=%g\n", norm.x(), norm.y(), norm.z(), norm.theta()/deg, norm.phi()/deg, mag);
484    if ( fabs(mag-1.)>kHalfTolerance ) printf("ERROR: not unitary.\n");
485    #endif
486
487    #ifdef SLANTWALLS
488    if ( isSlant && norm.z()==0. ) CorrectCylinderNormal(&norm);
489    #endif
490
491    return norm;
492}
493
494//__________________________________________________________________________________
495#ifdef SLANTWALLS
496void ShapeOfLens::CorrectCylinderNormal(G4ThreeVector* n) const {
497    n->setTheta( CLHEP::RandGauss::shoot(90.*deg, theSlantSigma) );
498    #ifdef DEBUG_SLANTWALLS
499        printf("%f\n",n->getTheta()/deg);
500    #endif
501}
502#endif
503
504//__________________________________________________________________________________
505G4double ShapeOfLens::GetCubicVolume(){
506  if(fCubicVolume == 0.)fCubicVolume=top_volume->GetCubicVolume();
507  return fCubicVolume;
508}
509
510//__________________________________________________________________________________
511G4double ShapeOfLens::GetSurfaceArea() {
512  if(fSurfaceArea == 0.)fSurfaceArea=top_volume->GetSurfaceArea();
513  return fSurfaceArea;
514}
515
516//__________________________________________________________________________________
517void ShapeOfLens::DescribeYourselfTo(G4VGraphicsScene& scene) const{
518    top_volume->DescribeYourselfTo(scene);
519}
520//__________________________________________________________________________________
521G4ThreeVector ShapeOfLens::GetPointOnSurface() const{
522    return top_volume->GetPointOnSurface();
523}
524
525//__________________________________________________________________________________
526G4bool ShapeOfLens::CalculateExtent(const EAxis pAxis,
527                           const G4VoxelLimits& pVoxelLimit,
528                           const G4AffineTransform& pTransform,
529                                 G4double& pmin, G4double& pmax) const {
530    return top_volume->CalculateExtent(pAxis,pVoxelLimit,pTransform,pmin,pmax);
531}
Note: See TracBrowser for help on using the repository browser.