source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/G4Detector/G4fresnellens/src/ArraySurface.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: 20.0 KB
Line 
1#include "ArraySurface.h"
2#include "FixedStack.h"
3
4#ifdef USE_ROOT
5#include <TROOT.h>
6#endif
7
8using namespace G4FresnelLens;
9
10int ArraySurface::fStepsLimit=100;
11#define NUM_TOLERANCE 1.e-15
12
13//______________________________________________________________________________
14ArraySurface::ArraySurface(int n, double* rho, double* z): SSurface(0., rho[n-1]), fN(n) {
15    fZ=new double [fN];
16    fRho=new double [fN];
17    memcpy(fZ, z, fN*sizeof(double));
18    memcpy(fRho, rho, fN*sizeof(double));
19
20    /*Nonfixed step*/// GetCurIndexPointer=&ArraySurface::GetCurIndexConstStep;
21    DetermineDistToSeg=&ArraySurface::DetermineDistToSegSmart;
22    //DetermineDistToSeg=&ArraySurface::DetermineDistToSegDirect;
23
24    //for (int i(0); i<fN; i++){
25        //printf("%i %g %g\n", i, fRho[i], fZ[i]);
26    //}
27
28    Configure();
29}
30
31//______________________________________________________________________________
32ArraySurface::~ArraySurface(){
33    delete [] fZ;
34    delete [] fRho;
35}
36
37//______________________________________________________________________________
38void ArraySurface::Configure(){
39    //TODO
40    fStep=fRho[1]-fRho[0];
41    fR1=fRho[0];
42
43    fZmin=fZmax=fZ[0];
44    for ( int i(1); i<fN; i++){
45        if ( fZ[i]<fZmin ) fZmin=fZ[i];
46        if ( fZ[i]>fZmax ) fZmax=fZ[i];
47    }
48    fZ1=fZ[0];
49    fZ2=fZ[fN-1];
50
51    double rho1=fRho[0];
52    double rho2=fRho[fN-1];
53    double z1=fZ[0];
54    double z2=fZ[fN-1];
55    int idxC=fN/2;
56    double rhoC=fRho[idxC];
57    double zC=z2-(rho2-rhoC)*(z2-z1)/(rho2-rho1);
58    double zCa=fZ[idxC];
59    fConvexity=zCa>zC?+1:-1;
60}
61
62//______________________________________________________________________________
63EInside ArraySurface::Inside(int idx, double rho, double z, FSIntersectedSegment *intSeg){
64    double zc1=fZ[idx];
65    double zc2=fZ[idx+1];
66    double rhoc1=fRho[idx];
67    double rhoc2=fRho[idx+1];
68    if ( intSeg ) FillSegment(idx, rho, intSeg);
69    return InsideCone(rho, z, rhoc1, zc1, rhoc2, zc2);
70}
71
72//______________________________________________________________________________
73double ArraySurface::DetermineDistToSegDirect(const G4ThreeVector& p, const G4ThreeVector& v, FSIntersectedSegment* intseg){
74    double dist_cone=kInfinity;
75
76    double perp=p.perp();
77    int idx=this->GetCurIndex(perp);
78    if( idx>=fN-1 ) idx=fN-2;
79    int found=-1;
80    double c=p.x()*v.x()+p.y()*v.y();
81    if ( c>=0. ) {
82        #ifdef DEBUG
83        if (theDebugSwitch) {
84            printf("go forward from %i to %i\n", idx, fN-1);
85        }
86        #endif
87        dist_cone=FindCone(p, v, idx, fN-1, +1, found);
88
89        #ifdef DEBUG
90        if (theDebugSwitch) {
91            printf("forward dist=%f and idx=%i\n", dist_cone, found);
92        }
93        #endif
94    }
95    else {
96        double diff=p.perp2()-c*c/v.perp2();
97        int min_i=this->GetCurIndex(diff>0.?sqrt(diff):0.);
98
99        #ifdef DEBUG
100        if (theDebugSwitch) {
101            printf("current i1=%i min_i=%i perp=%f\n", idx, min_i, perp);
102        }
103        #endif
104
105        if ( min_i<0) min_i=0;
106        if ( min_i<idx) {
107            #ifdef DEBUG
108            if (theDebugSwitch) {
109                printf("go backward from %i to %i\n",idx,0);
110            }
111            #endif
112            dist_cone=FindCone(p, v, idx, min_i-1, -1, found);
113            #ifdef DEBUG
114            if (theDebugSwitch) {
115                printf("backward dist=%f and idx=%i\n",dist_cone,found);
116            }
117            #endif
118        }
119        if ( dist_cone==kInfinity ) {
120            #ifdef DEBUG
121            if (theDebugSwitch) {
122                printf("go forward from %i to %i\n",idx,fN-1);
123            }
124            #endif
125            dist_cone=FindCone(p, v, idx, fN-1, +1, found);
126            #ifdef DEBUG
127            if (theDebugSwitch) {
128                printf("forward dist=%f and idx=%i\n",dist_cone,found);
129            }
130            #endif
131        }
132    }
133
134    if ( intseg && found>-1 ){
135        FillSegment(found, -1, intseg);
136    }
137    return dist_cone;
138}
139
140//______________________________________________________________________________
141double ArraySurface::DistanceToSurfSafe(const G4ThreeVector& p){
142    double perp2=p.perp2();
143    if ( perp2<=fR*fR && perp2>=fR1*fR1 ) {
144        if ( p.z()>fZmax ) return p.z()-fZmax;
145        if ( p.z()<fZmin ) return fZmin-p.z();
146        return 0.;
147    }
148    double perp=sqrt(perp2);
149    double rc=perp>fR?fR:fR1;
150    if ( p.z()>fZmax ) return Hyp(perp, p.z(), rc, fZmax);
151    if ( p.z()<fZmin ) return Hyp(perp, p.z(), rc, fZmin);
152    return fabs(rc-perp);
153}
154
155//__________________________________________________________________________________
156double ArraySurface::FindCone(const G4ThreeVector& p, const G4ThreeVector&  v, int point1, int point2, int dir, int& idx){
157    double dist;
158    #ifdef DEBUG
159        if (theDebugSwitch) {
160            printf("pos %e %e %e\n", p.x(), p.y(), p.z());
161            printf("dir %e %e %e\n", v.x(), v.y(), v.z());
162        }
163    #endif
164    if ( (point1-point2)*dir>0 ) Abort_("Bad iteration direction. Shouldn't happen.");
165    for (int i(point1); i!=point2; i+=dir){
166        dist=ConeIntersection(p, v, fRho[i], fZ[i], fRho[i+1], fZ[i+1]);
167        if ( dist<kInfinity ){
168            idx=i;
169            #ifdef DEBUG
170                if (theDebugSwitch) {
171                    info_<<"Direct alg: yes "<<i-point1<<endline_;
172                    info_<<"from/to/found "<<point1<<"  "<<point2<<"  "<<i<<endline_;
173                }
174            #endif
175            return dist;
176        }
177    }
178
179    #ifdef DEBUG
180        if (theDebugSwitch) {     info_<<"Direct alg: noo "<<point2-point1<<endline_;
181            info_<<"from/to "<<point1<<"  "<<point2<<endline_;
182        }
183    #endif
184    idx=-1;
185    return kInfinity;
186}
187
188//__________________________________________________________________________________
189double ArraySurface::FindConeConcave(double aperp, const G4ThreeVector& p, const G4ThreeVector&  v, int aidx, FSIntersectedSegment* intseg){
190    G4ThreeVector pos(p);
191    double perp=aperp;
192    int idx=-1;
193    double dist(kInfinity);
194    int currentI=aidx;
195    bool out=VectorGoesOut(p,v);
196
197    fixedISStack stackSteps;
198    //stackSteps.clear();
199    if ( out ){
200        if ( currentI>0 ) stackSteps.push(0, currentI);   //left
201        if ( (fN-1-currentI)>1 ) stackSteps.push(currentI+1, fN-1); //right
202    }
203    else {
204        if ( (fN-1-currentI)>1 ) stackSteps.push(currentI+1, fN-1); //right
205        if ( currentI>0 ) stackSteps.push(0, currentI);   //left
206    }
207    stackSteps.push(currentI, currentI+1);                            //center
208
209    #ifdef DEBUG
210    if (theDebugSwitch) {
211        printf("Starting concave method for\n");
212        info_<<"currentI "<<currentI<<endline_;
213    }
214    #endif
215
216    fcIterationStep* itStep;
217    for ( int i(0); i<fStepsLimit; i++ ){
218        itStep = stackSteps.pop();
219        if ( !itStep ) {
220            #ifdef DEBUG
221            if (theDebugSwitch) {
222                info_<<"Concave alg: no "<<i<<endline_;
223            }
224            #endif
225            return kInfinity;
226        }
227
228        #ifdef DEBUG
229        if (theDebugSwitch) {
230            info_<<i<<endline_;
231            info_<<Form(" min %i max %i c %i", 0, fN-1, currentI)<<endline_;
232            info_<<Form("tt l %f l1 %f ", fRho[0], fRho[fN-1])<<endline_;
233            info_<<Form("coord %s (%e %e) (%e %e) idx(%i, %i) %e %e %e %e", (out?"out":"inn"),
234                         pos.x(), pos.y(), pos.perp(), pos.z(),
235                         itStep->leftI, itStep->rightI,
236                         fRho[itStep->leftI], fZ[itStep->leftI], fRho[itStep->rightI], fZ[itStep->rightI])<<endline_;
237        }
238        #endif
239        dist=ConeIntersection( p, v, fRho[itStep->leftI], fZ[itStep->leftI], fRho[itStep->rightI], fZ[itStep->rightI] );
240
241        #ifdef DEBUG
242        if ( theDebugSwitch ) {
243            info_<<Form("dist=%g", dist)<<endline_;
244            info_<<Form("dir %e %e %e", v.x(), v.y(), v.z())<<endline_;
245        }
246        #endif
247        if ( !(dist<kInfinity) ) continue;
248
249        //find intersection position
250        pos=p+v*dist;
251        perp=pos.perp();
252        int backi=currentI;
253        currentI=this->GetCurIndex(perp);
254
255        if ( (itStep->rightI-itStep->leftI)==1 ){
256            idx=currentI;
257            #ifdef DEBUG
258            if (theDebugSwitch) {
259                info_<<"Concave alg: yes "<<i<<endline_;
260            }
261            #endif
262            if (intseg) FillSegment(idx, perp, intseg);
263            return dist;
264        }
265
266        // calculate new steps
267        int li(itStep->leftI), ri(itStep->rightI);
268        if ( currentI==backi ) {
269            if ( currentI==li-1 ) {
270                currentI=li;
271                if (li) li--;
272            }
273            else if ( currentI==ri+1 ) {
274                currentI=ri;
275                if ( ri<fN-1 ) ri++;
276            }
277            #ifdef DEBUG
278            else {
279                printf("strange ci\n");
280            }
281            #endif
282        }
283        #ifdef DEBUG
284        if ( theDebugSwitch )
285          printf("add single seg at %i (prev %i)\n", currentI, backi);
286        #endif
287
288        if ( VectorGoesOut(pos,v) ){
289          if ( (currentI-li)>0   )
290            stackSteps.push(li, currentI);   //left
291          if ( (ri-currentI)>1 )
292            stackSteps.push(currentI+1, ri); //right
293        }
294        else {
295          if ( (ri-currentI)>1 )
296            stackSteps.push(currentI+1, ri); //right
297          if ( (currentI-li)>0   )
298            stackSteps.push(li, currentI);   //left
299        }
300        stackSteps.push(currentI, currentI+1);
301    }
302    error_<<"FindConeConcave returned no result after "<<fStepsLimit<<" steps"<<endline_;
303    #ifdef DEBUG
304        throw 1;
305    #endif
306    return kInfinity;
307}
308
309//__________________________________________________________________________________
310double ArraySurface::FindConeConvex(double aperp, const G4ThreeVector& p, const G4ThreeVector&  v, int aidx, FSIntersectedSegment* intseg){
311    G4ThreeVector pos(p);
312    int idx=-1;
313    int currentI=aidx;
314    double perp2(p.perp2()), perp(aperp);
315    double dist;
316    double minRho=fRho[0];
317    double maxRho=fRho[fN-1];
318    double minZ=fZmin;
319    double maxZ=fZmax;
320
321    bool inlimits=false;
322
323    #ifdef DEBUG
324        if (theDebugSwitch) {
325            printf("Starting convex method\n");
326            printf("initins perp=%e z=%e r1=%e z1=%e r2=%e z2=%e %i\n",
327            perp, pos.z(),
328            fRho[currentI], fZ[currentI],
329            fRho[currentI+1], fZ[currentI+1], currentI);
330        }
331    #endif
332
333    double mindist=0.;
334    for (int i(0); i<=fStepsLimit; i++){
335        #ifdef DEBUG
336        if (theDebugSwitch) {
337            info_<<"Iteration "<<i<<endline_;
338            info_<<"i cur/min/max "<<currentI<<"  "<<0<<"   "<<fN-1<<endline_;
339            info_<<"rho cur/min/max "<<fRho[currentI]<<"  "<<fRho[0]<<"   "<<fRho[fN-1]<<endline_;
340        }
341        #endif
342
343        dist=ConeIntersection1( p, v, fRho[currentI], fZ[currentI], fRho[currentI+1], fZ[currentI+1], inlimits, mindist);
344        #ifdef DEBUG
345        if (theDebugSwitch) {
346            double dist1=ConeIntersection( p, v, fRho[currentI], fZ[currentI], fRho[currentI+1], fZ[currentI+1] );
347            bool inlimits1;
348            double dist2=ConeIntersection1( p, v, fRho[currentI], fZ[currentI], fRho[currentI+1], fZ[currentI+1], inlimits1, 0 );
349            info_<<Form("dist_cone=%e (raw %e), (wo mindist %e) in limits=%i (%i)", dist, dist1, dist2, inlimits, inlimits1)<<endline_;
350            printf("coord (%e %e) (%e %e) %e %e %e %e\n",
351                    pos.x(), pos.y(), pos.perp(), pos.z(),
352                    fRho[currentI], fZ[currentI], fRho[currentI+1], fZ[currentI+1]);
353            printf("dir %g %g\n", v.theta()/deg, v.phi()/deg);
354        }
355        #endif
356
357        if ( dist==kInfinity ) {
358            #ifdef DEBUG
359            if (theDebugSwitch) {
360                info_<<"Conv alg: noo "<<i<<endline_;
361            }
362            #endif
363            return kInfinity;
364        }
365        else if( inlimits ) {
366            idx=currentI;
367            #ifdef DEBUG
368            if (theDebugSwitch) {
369                info_<<"Conv alg: yes "<<i<<endline_;
370            }
371            #endif
372            if (intseg) FillSegment(idx, -1., intseg);
373            return dist;
374        }
375        mindist=dist*0.9; //FIXME: not optimized
376        pos=p+v*dist;
377        perp2=pos.perp2();
378
379        #ifdef DEBUG
380         if (theDebugSwitch) {
381            printf("perp=%e z=%e r1=%e z1=%e r2=%e z2=%e i1=%i i2=%i\n", sqrt(perp2), pos.z(), minRho, minZ, maxRho, maxZ, 0, fN-1);
382         }
383        #endif
384        if ( !InLimits2(perp2, pos.z(), minRho, minZ, maxRho, maxZ) ) {
385            #ifdef DEBUG
386            if (theDebugSwitch) {
387                info_<<Form("Conv alg: noo not in limits %i: perp=%e z=%e r1=%e z1=%e r2=%e z2=%e", i, sqrt(perp2), pos.z(), minRho, minZ, maxRho, maxZ)<<endline_;
388            }
389            #endif
390            return kInfinity;
391        }
392        perp=sqrt(perp2);
393        currentI=this->GetCurIndex(perp);
394        if ( currentI==fN-1 ){
395            #ifdef DEBUG
396            if (theDebugSwitch) {
397                info_<<"Conv alg: noo "<<i<<endline_;
398             }
399            #endif
400            return kInfinity;
401        }
402
403        #ifdef DEBUG
404        if (theDebugSwitch) {
405            debug_<<currentI<<endline_;
406        }
407        #endif
408    }
409    error_<<"FindConeConvex returned no result after "<<fStepsLimit<<" steps."<<endline_;
410    #ifdef DEBUG
411        throw 1;
412    #endif
413    return kInfinity;
414}
415
416//__________________________________________________________________________________
417double ArraySurface::ConeIntersection(const G4ThreeVector& p, const G4ThreeVector& v, double r1, double z1, double r2, double z2){
418    // adopted from Ltrace Lrt_bsearch
419    // Int_t Lcheck_intersect_cone(Double_t r1, Double_t z1, Double_t r2, Double_t z2,
420    //                Double_t phot[8], Double_t *dis)
421//    printf("(%10.15f, %10.15f, %10.15f), (%10.15f, %10.15f, %10.15f)\n", p.x(), p.y(), p.z(), v.x(), v.y(), v.z());
422//    printf("(%10.15f, %10.15f), (%10.15f, %10.15f)\n", r1, r2, z1, z2);
423
424    double d1, d2;
425    int nroot;
426    if ( z1==z2 ) {
427        // the quadratic equation for this case fails due to numeric precision
428        double vz=v.z();
429        if ( vz ) {
430            d1=d2=(z1-p.z())/v.z();
431            nroot=1;
432        }
433        else return kInfinity;
434    }
435    else {
436        double a00 = (z2-z1)/(r2-r1);
437        double z0  = z1-a00*r1;
438        double a0  = a00*a00;
439
440        double dz = p.z()-z0;
441        double a = a0*v.perp2() - v.z()*v.z();
442        double b = 2.0*(a0*v.x()*p.x() + a0*v.y()*p.y() - v.z()*dz);
443        double c = a0*p.perp2() - dz*dz;
444        nroot=SolveQuadratic(a, b, c, d1, d2);
445    }
446
447    //    if ( r1>=937.00 && r1<=938.0) printf("nroot=%i d1=%e d2=%e\n",nroot,d1,d2);
448    //there can be 1 intersection
449    if( nroot==0 || (nroot==1&&d1==0.0) ) return kInfinity;
450    // d2 is more than d1 as defined in SolveQuadratic
451    if(d2<0.0)  return kInfinity;
452
453    G4ThreeVector p1;
454    double z;
455    if ( d1>=0.0 ){
456        p1=p+v*d1;
457        z=p1.z();
458        /* check if the solution is in the specified range */
459//        printf("dd %e (%10.15f, %10.15f, %10.15f)\n", d1, p1.x(), p1.y(), p1.z() );
460        if( (z-z1)*(z-z2)<=NUM_TOLERANCE ) {
461            double perp2=p1.perp2();
462            if ( perp2>=r1*r1 && perp2<=r2*r2 ) return d1;
463        }
464    }
465
466    if ( d1!=d2 ) {
467        p1=p+v*d2;
468        z=p1.z();
469        /* check if the solution is in the specified range */
470//        printf("dd %e (%10.15f, %10.15f, %10.15f)\n", d2, p1.x(), p1.y(), p1.z() );
471        if( (z-z1)*(z-z2)<=NUM_TOLERANCE ) {
472            double perp2=p1.perp2();
473            if ( perp2>=r1*r1 && perp2<=r2*r2 ) return d2;
474        }
475    }
476    return kInfinity;
477}
478
479//__________________________________________________________________________________
480double ArraySurface::ConeIntersection1(const G4ThreeVector& p, const G4ThreeVector& v,
481                                    double r1, double z1, double r2, double z2,
482                                    bool& inlimits, double mindist){
483    // adopted from Ltrace Lrt_bsearch
484    // Int_t Lcheck_intersect_cone(Double_t r1, Double_t z1, Double_t r2, Double_t z2,
485    //                Double_t phot[8], Double_t *dis)
486//    printf("(%10.15f, %10.15f, %10.15f), (%10.15f, %10.15f, %10.15f)\n", p.x(), p.y(), p.z(), v.x(), v.y(), v.z());
487//    printf("(%10.15f, %10.15f), (%10.15f, %10.15f)\n", r1, r2, z1, z2);
488    inlimits=false;
489    double d1, d2, z0;
490    int nroot;
491    double cone_higher; // help to chose right cone part and reject intersections with the reflected cone
492    if ( z1==z2 ) {
493        // the quadratic equation for this case fails due to numeric precision
494        double vz=v.z();
495        if ( vz ) {
496            d1=d2=(z1-p.z())/v.z();
497            nroot=1;
498        }
499        else return kInfinity;
500
501        z0=z2;
502        cone_higher=0.;
503    }
504    else {
505        double a00 = (z2-z1)/(r2-r1);
506        z0  = z1-a00*r1;
507        cone_higher=(z2-z0);
508        double a0  = a00*a00;
509
510        double dz = p.z()-z0;
511        double a = a0*v.perp2() - v.z()*v.z();
512        double b = 2.0*(a0*v.x()*p.x() + a0*v.y()*p.y() - v.z()*dz);
513        double c = a0*p.perp2() - dz*dz;
514
515        nroot=SolveQuadratic(a, b, c, d1, d2);
516    }
517
518    //    if ( r1>=937.00 && r1<=938.0) printf("nroot=%i d1=%e d2=%e\n",nroot,d1,d2);
519    //there can be 1 intersection
520    if( nroot==0 || (nroot==1&&d1==0.0) ) return kInfinity;
521    // d2 is more than d1 as defined in SolveQuadratic
522    if(d2<0.0) return kInfinity;
523    #ifdef CHECK_ConeIntersection1
524    printf("dd000 %.10f %.10f\n", d1, d2);
525    #endif
526
527    G4ThreeVector p1;
528    double z;
529    double minpositive=kInfinity;
530    if ( d1>=mindist ){
531        p1=p+v*d1;
532        z=p1.z();
533        #ifdef CHECK_ConeIntersection1
534        printf("dd %e (%10.15f, %10.15f), (pos %10.15f, %10.15f, %s), (cone %10.15f, %10.15f, %s)\n", d1, p1.perp(), p1.z(), z, z0, ((z>=z0)?"aga":"nee"), z1, z0, (cone_higher>=0?"aga":"nee") );
535        #endif
536        // check that solution is on the same halfspace
537        if ( (z-z0)*cone_higher>=0. ) {
538            minpositive=d1;
539            if( (z-z1)*(z-z2)<NUM_TOLERANCE ) {
540                double perp2=p1.perp2();
541                if ( perp2>=r1*r1 && perp2<=r2*r2 ) {
542                    inlimits=true;
543                    return d1;
544                }
545            }
546        }
547        //else minpositive=d2;
548    }
549
550    if ( d1!=d2 && d2>=mindist ){
551        p1=p+v*d2;
552        z=p1.z();
553        /* check if the solution is in the specified range */
554        #ifdef CHECK_ConeIntersection1
555        printf("dd %e (%10.15f, %10.15f), (pos %10.15f, %10.15f, %s), (cone %10.15f, %10.15f, %s)\n", d2, p1.perp(), p1.z(), z, z0, ((z>=z0)?"aga":"nee"), z1, z0, (cone_higher?"aga":"nee") );
556        #endif
557        if ( (z-z0)*cone_higher>=0. ) {
558            if( (z-z1)*(z-z2)<NUM_TOLERANCE ) {
559                double perp2=p1.perp2();
560                if ( perp2>=r1*r1 && perp2<=r2*r2 ) {
561                    inlimits=true;
562                    return d2;
563                }
564            }
565
566            if ( d2<minpositive ) minpositive=d2;
567        }
568        // else minpositive=kInfinity; //FIXME
569    }
570
571    inlimits=false;
572    return minpositive;
573}
574
575//__________________________________________________________________________________
576EInside ArraySurface::InsideCone(double rho, double z, double rho1, double z1, double rho2, double z2){
577    double z0=(z2-z1)*(rho-rho1)/(rho2-rho1)+z1;
578
579    // find a distance from the point to the cylinder surface using the triangle area
580    double s=(rho-rho1)*(z2-z1)-(rho2-rho1)*(z-z1);
581    double drho=rho2-rho1;
582    double dz=z2-z1;
583    double h2=s/sqrt(drho*drho+dz*dz);
584    if ( h2>kHalfTolerance ){
585        if ( z>z0 ) return kOutside;
586        else return kInside;
587    }
588    return kSurface;
589}
590
591//__________________________________________________________________________________
592bool ArraySurface::operator==(SSurface& s){
593    ArraySurface* to1=static_cast<ArraySurface*>(&s);
594    if (!to1) {
595        return false;
596    }
597    ArraySurface& to=*to1;
598    if( fN!=to.fN ||
599        fStep!=to.fStep ||
600        fConvexity!=to.fConvexity ||
601        fZmin!=to.fZmin ||
602        fZmax!=to.fZmax ){
603        //print();
604        //to.print();
605        return false;
606    }
607    for (int i(0); i<fN; i++){
608        if ( fZ[i]==to.fZ[i] &&
609             fRho[i]==to.fRho[i] ) continue;
610        //printf("%g!=%g, %g!=%g\n", fZ[i], to.fZ[i], fRho[i], to.fRho[i]);
611        return false;
612    }
613    return true;
614}
Note: See TracBrowser for help on using the repository browser.