1 | #include "ArraySurface.h" |
---|
2 | #include "FixedStack.h" |
---|
3 | |
---|
4 | #ifdef USE_ROOT |
---|
5 | #include <TROOT.h> |
---|
6 | #endif |
---|
7 | |
---|
8 | using namespace G4FresnelLens; |
---|
9 | |
---|
10 | int ArraySurface::fStepsLimit=100; |
---|
11 | #define NUM_TOLERANCE 1.e-15 |
---|
12 | |
---|
13 | //______________________________________________________________________________ |
---|
14 | ArraySurface::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 | //______________________________________________________________________________ |
---|
32 | ArraySurface::~ArraySurface(){ |
---|
33 | delete [] fZ; |
---|
34 | delete [] fRho; |
---|
35 | } |
---|
36 | |
---|
37 | //______________________________________________________________________________ |
---|
38 | void 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 | //______________________________________________________________________________ |
---|
63 | EInside 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 | //______________________________________________________________________________ |
---|
73 | double 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 | //______________________________________________________________________________ |
---|
141 | double 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 | //__________________________________________________________________________________ |
---|
156 | double 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 | //__________________________________________________________________________________ |
---|
189 | double 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 | //__________________________________________________________________________________ |
---|
310 | double 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 | //__________________________________________________________________________________ |
---|
417 | double 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 | //__________________________________________________________________________________ |
---|
480 | double 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 | //__________________________________________________________________________________ |
---|
576 | EInside 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 | //__________________________________________________________________________________ |
---|
592 | bool 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 | } |
---|