source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/G4Detector/G4fresnellens/src/FresnelSurface.cc @ 114

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

actual version of ESAF at CCin2p3
File size: 11.1 KB
Line 
1#include "FresnelSurface.h"
2#include "FresnelOutput_def.h"
3#include "LensUtils.h"
4#include "CylindricSurface.h"
5
6#include <CLHEP/Matrix/Matrix.h>
7using CLHEP::HepMatrix;
8
9#include <G4GeometryTolerance.hh>
10#include <globals.hh>
11
12#include <float.h>
13#include <stdarg.h>
14
15using namespace G4FresnelLens;
16using std::swap;
17using std::min;
18using std::max;
19using std::vector;
20
21//__________________________________________________________________________________
22FresnelSurface::FresnelSurface(const std::vector<SSurface*>& surf,
23                               ParametricSurface* downsurf,
24                               ParametricSurface* upsurf):
25SSurface(),
26fUpLimit(upsurf),
27fDownLimit(downsurf),
28theToothEdges(0)
29{
30    FillArray(surf);
31}
32
33//__________________________________________________________________________________
34FresnelSurface::~FresnelSurface(){
35    if (theToothEdges) delete [] theToothEdges;
36    if ( fUpLimit ) delete fUpLimit;
37    if ( fDownLimit ) delete fDownLimit;
38
39    for (unsigned i(0); i<fSurfaces.size(); i++){
40        delete fSurfaces[i];
41    }
42    fSurfaces.clear();
43
44    for (unsigned i(0); i<fSurfacesIns.size(); i++){
45        if (i%2) delete fSurfacesIns[i];
46    }
47    fSurfacesIns.clear();
48}
49
50//__________________________________________________________________________________
51void FresnelSurface::FillArray(const std::vector<SSurface*>& surf){
52    #ifdef DEBUG
53      theDebugSwitch=false;
54    #endif
55    theLastPerp=-1.;
56    theLastInside=(EInside)0;
57    if (theToothEdges) delete [] theToothEdges;
58
59    fR=surf.back()->GetR2();
60    fUpLimit->SetR(fR);
61    fUpLimit->ShiftSurfaceTo(surf, false);
62    fDownLimit->SetR(fR);
63    fDownLimit->ShiftSurfaceTo(surf, true);
64
65    fCylinder.SetR(fR);
66    fCylinder.SetZ(fDownLimit->GetMinZ(), fUpLimit->GetMaxZ());
67
68    theToothN=surf.size();
69    theToothEdges = new double[theToothN+1];
70    theToothEdges[0]=0.;
71
72    for (unsigned i(0); i<surf.size()-1; i++){
73        double r=surf[i]->GetR();
74        theToothEdges[i+1]=r;
75        fSurfaces.push_back(surf[i]);
76        fSurfaces.push_back(
77            new CylindricSurface(surf[i]->GetZ2(),
78                                 surf[i+1]->GetZ1(), r));
79        ((CylindricSurface*)fSurfaces.back())->SetZ(fDownLimit->PointOnSurf(r), fUpLimit->PointOnSurf(r));
80        fSurfacesIns.push_back(surf[i]);
81        fSurfacesIns.push_back(
82            new CylindricSurface(surf[i]->GetZ2(),
83                                 surf[i+1]->GetZ1(), r));
84    }
85    fSurfaces.push_back(surf[surf.size()-1]);
86    fSurfacesIns.push_back(surf[surf.size()-1]);
87    theToothEdges[surf.size()]=surf.back()->GetR();
88
89//     #define PRINT
90    #ifdef PRINT
91    fUpLimit->print();
92    fDownLimit->print();
93    for (unsigned int i(0); i<fSurfaces.size(); i++){
94        printf("%5u  ", i); fSurfaces[i]->print();
95        if (i%2) {printf("%5u* ", i); fSurfacesIns[i]->print();}
96    }
97    #endif
98}
99
100//__________________________________________________________________________________
101#ifdef USE_ROOT
102#include <TGraph.h>
103#include <TCanvas.h>
104#include <TString.h>
105#include <TFile.h>
106void FresnelSurface::draw(const char* file, const char* canvasfile, double shift, const char* prefix, const char* option){
107    TString file1=file;
108    TFile* f = new TFile(file1, option);
109    TGraph surf, ulim, dlim;
110    ulim.SetLineColor(kGreen);
111    dlim.SetLineColor(kRed);
112
113    vector<SSurface*>::iterator it=fSurfaces.begin();
114    for ( ; it<fSurfaces.end(); it+=2 ){
115        int i(0);
116        double rho,z;
117        while((*it)->Point(i, 0., rho, z)) {
118            surf.SetPoint(surf.GetN(), rho, z+shift);
119            i++;
120        }
121    }
122
123    double step=0.01, rho, z;
124    int i(0);
125    while( fUpLimit->Point(i, step, rho, z) ){
126        ulim.SetPoint(surf.GetN(), rho, z+shift);
127        i++;
128    }
129    i=0;
130    while( fDownLimit->Point(i, step, rho, z) ){
131        ulim.SetPoint(surf.GetN(), rho, z+shift);
132        i++;
133    }
134    surf.Write( TString(prefix)+"surf" );
135    ulim.Write( TString(prefix)+"up" );
136    dlim.Write( TString(prefix)+"down" );
137    f->Close();
138
139    if (canvasfile) {
140        TCanvas c("surface", "surface");
141        surf.Draw("AL");
142        ulim.Draw("L");
143        dlim.Draw("L");
144        c.Print(canvasfile);
145    }
146}
147#endif
148
149//__________________________________________________________________________________
150EInside FresnelSurface::Inside(double rho, double z, FSIntersectedSegment *intSeg){
151    #ifdef DEBUG
152      if ( theDebugSwitch ){
153          printf("FresnelSurface::Inside (%.10f, %.10f)\n", rho, z);
154      }
155    #endif
156    if ( fabs(theLastPerp-rho)<kHalfTolerance && fabs(theLastZ-rho)<kHalfTolerance ) {
157         #ifdef DEBUG
158           if ( theDebugSwitch ){
159               printf( "Return last saved state (%.10f, %.10f): %s\n", rho, z,
160                       theLastInside==kSurface?"surface":(theLastInside==kInside?"inside":"outside"));
161           }
162         #endif
163         if ( intSeg && theLastInside==kSurface ) *intSeg=theLastIntSeg;
164         return theLastInside;
165    }
166    theLastPerp=rho;
167    theLastZ=z;
168
169    #ifdef DEBUG
170      if ( theDebugSwitch ){
171          printf("Check top/bottom coordinates: b(%.10f), z(%.10f), t(%.10f)\n", this->GetBottomZ(rho), z, this->GetUpperZ(rho) );
172      }
173    #endif
174    if ( z<GetBottomZ(rho)) return theLastInside=kInside;
175    if ( z>GetUpperZ(rho) ) return theLastInside=kOutside;
176
177    #ifdef DEBUG
178      if ( theDebugSwitch ){
179          printf("Find tooth\n");
180      }
181    #endif
182    int n_tooth = BinarySearch(theToothEdges,rho,0,theToothN);
183    if( n_tooth<0 ) {
184        printf("Warning. Tooth not found. Shouldn't happen (rho=%g of [%g, %g]).",
185               rho, theToothEdges[0], theToothEdges[theToothN]);
186        return theLastInside=kOutside;
187    }
188
189    vector<SSurface*>::const_iterator it=fSurfacesIns.begin()+2*n_tooth;
190    if ( (n_tooth                 && (*(it-1))->Inside(rho, z, intSeg)==kSurface) || //left
191         (it+1<fSurfacesIns.end() && (*(it+1))->Inside(rho, z, intSeg)==kSurface)    //right
192       ) {
193           if (intSeg) theLastIntSeg=*intSeg;
194           return theLastInside=kSurface;
195       }
196    theLastInside=(*it)->Inside(rho, z, intSeg);                              //surface
197    if ( intSeg && theLastInside==kSurface ) theLastIntSeg=*intSeg;
198    return theLastInside;
199}
200
201//__________________________________________________________________________________
202double FresnelSurface::FindNearestIntersection(const G4ThreeVector& p, const G4ThreeVector& v, FSIntersectedSegment *intSeg){
203    //
204    // find nearest intersection with fresnel surface
205    //
206    #ifdef DEBUG
207      if (theDebugSwitch) {
208          info_<<"FindeNearestIntersection started"<<endline_;
209      }
210    #endif
211
212    // initialize tracking variables
213    G4ThreeVector pos=p;
214    double perp=p.perp();
215    double shift=0.;
216
217    {// approach to the limiting surface
218        #ifdef DEBUG
219          if (theDebugSwitch) {
220              info_<<"Approaching"<<endline_;
221          }
222        #endif
223        double dist(kInfinity);
224        double z_down=this->GetBottomZ(perp);
225        if ( pos.z()<z_down ) dist=fDownLimit->DistanceToSurf(pos, v);
226        else {
227            double z_up=this->GetUpperZ(perp);
228            if ( pos.z()>z_up ) {
229                dist=fUpLimit->DistanceToSurf(pos, v);
230            }
231        }
232        #ifdef DEBUG
233        if (theDebugSwitch) {
234            printf("ddd %e z %e zdown %e\n", dist, pos.z(), z_down);
235        }
236        #endif
237        if ( dist==kInfinity ) {
238            if ( perp>fR ){
239                dist=fCylinder.DistanceToSurf(pos, v);
240                if ( dist==kInfinity ) return kInfinity;
241                pos+=v*dist;
242                shift+=dist;
243            }
244        }
245        else {
246            pos+=v*dist;
247            shift+=dist;
248        }
249        perp=pos.perp();
250        #ifdef DEBUG
251        if (theDebugSwitch) {
252            printf("perp=%f shift=%f z=%f\n", perp, shift, pos.z());
253        }
254        #endif
255    }// approach to the limiting surface
256
257    int theCurrentToothIndex=BinarySearch(theToothEdges, perp, 0, theToothN);
258    if ( theCurrentToothIndex<0 ) return kInfinity;
259    vector<SSurface*>::const_iterator it=fSurfaces.begin()+theCurrentToothIndex*2;
260    #ifdef DEBUG
261    if (theDebugSwitch) {
262        printf("i=%i rhotooth=%f\n", theCurrentToothIndex, theToothEdges[theCurrentToothIndex]);
263    }
264    #endif
265
266    while ( true ){
267        if ( intSeg ) intSeg->Reset();
268        int direction;
269        {// determine next tooth index
270            double pv=pos.x()*v.x()+pos.y()*v.y();
271            double c2(0.), c2a(0.);
272            //printf("%g pv\n", pv);
273            if ( pv>=0. ) {
274                direction=+1;
275            }
276            else {
277                double ip2=1./pos.perp2();
278                double tp=(*it)->GetR1();
279                //printf("tp %g\n", tp);
280                if ( tp==0. ) direction=+1.;
281                else {
282                    double r_over_rho=tp*tp*ip2;
283                    if ( r_over_rho>1. ) r_over_rho=1.;
284                    c2a=r_over_rho-1.;
285                    c2=(pv>=0.?1.:-1.)*pv*pv*ip2/v.perp2();
286                    //printf("c2 %g c2a %g\n", c2, c2a);
287                    if ( c2<-1. ) c2=-1.;
288                    direction = c2>=c2a?+1:-1;
289                }
290            }
291        }// determine next tooth index
292        vector<SSurface*>::const_iterator next=it+direction;
293
294        FSIntersectedSegment cylseg;
295        double dist_cyl = next<fSurfaces.end() ? (*next)->DistanceToSurf(pos, v, &cylseg) : kInfinity;
296        double dist_cone=(*it)->DistanceToSurf(pos, v, intSeg);
297        if ( dist_cone<dist_cyl ) {
298            if ( ( (dist_cyl-dist_cone)<=kTolerance ) ) {
299                // skip edge with the thickness less than kTolerance,
300                it+=2*direction;
301                if ( it>=fSurfaces.end() ) return kInfinity;
302                shift+=dist_cyl; pos+=v*dist_cyl; perp=pos.perp();
303                continue;
304            }
305            return shift+dist_cone;
306        }
307        else if ( dist_cyl==kInfinity ) return kInfinity;
308
309        {//shift to cylinder
310            int step=theCurrentToothIndex*2+direction; 
311            if ( step<0 || step>=(int)fSurfacesIns.size() ) {
312                printf("Trying to step to the unexistent cylinder:\n");
313                printf("idx=%i, dir=%i, step=%i, size=%i\n", theCurrentToothIndex, direction, step, (int)fSurfacesIns.size());
314                return kInfinity;
315            }
316
317            if ( dist_cyl<kHalfTolerance ) dist_cyl=kHalfTolerance;
318            shift+=dist_cyl; pos+=v*dist_cyl; perp=pos.perp();
319
320            vector<SSurface*>::const_iterator it1=fSurfacesIns.begin()+step;
321            if ( (pos.z()-(*it1)->GetZ1())*(pos.z()-(*it1)->GetZ2())<0 ){
322                if ( (*(it1+direction))->TestEdgeTolerance(direction, pos, v)>kTolerance ){
323                    if ( intSeg ) *intSeg=cylseg;
324                    return shift;
325                }
326                shift+=kTolerance; pos+=v*kTolerance; perp=pos.perp();
327            }
328        }//shift to cylinder
329        it+=2*direction;
330        if ( it>=fSurfaces.end() ) return kInfinity;
331    }
332}
333
334//______________________________________________________________________________
335double* FresnelSurface::CopyTeeth(int& n){
336    double* teeth=new double [theToothN+1];
337    n=theToothN+1;
338    memcpy(teeth, theToothEdges, (theToothN+1)*sizeof(double));
339    return teeth;
340}
Note: See TracBrowser for help on using the repository browser.