source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/optics/src/DetectorPhotonTransporter.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.4 KB
Line 
1// $Id: DetectorPhotonTransporter.cc 2918 2011-06-10 22:22:31Z mabl $
2// Author: D.Demarco, M.Pallavicini
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: DetectorPhotontransporter                                            *
8 *  Package: Optics                                                          *
9 *  Coordinator: Alessandro.Thea                                             *
10 *                                                                           *
11 *****************************************************************************/
12
13//______________________________________________________________________________
14//
15//  DetectorPhotonTransporter
16//  =========================
17//
18//  Abstract base class providing the interface of a generic photon
19//  transporter. Every element in the detector that is able to interact
20//  with photons is a DetectorPhotonTransporter.
21//  The DetectorPhotonTransporters are represented as cylinders with position
22//  and height (fDZdown+fDZup).
23//
24
25#include "DetectorPhotonTransporter.hh"
26#include "Etypes.hh"
27#include <TMath.h>
28ClassImp(DetectorPhotonTransporter)
29
30using namespace TMath;
31
32//______________________________________________________________________________
33DetectorPhotonTransporter::DetectorPhotonTransporter(): EsafConfigurable(),
34    fPos(0,0,0), fDZdown(0), fDZup(0), fR(0) {
35    //
36    // Constructor
37    //
38
39//    fLastCylInt.first=0;
40//    fLastCylInt.second=new Double_t[2];
41}
42
43//______________________________________________________________________________
44DetectorPhotonTransporter::~DetectorPhotonTransporter(){
45    //
46    // Destructor
47    //
48
49//    delete [] fLastCylInt.second;
50}
51
52
53//______________________________________________________________________________
54Bool_t DetectorPhotonTransporter::IsInside( Photon *p ) const {
55    //
56    // True if the photon is inside the boundaries of the transporter
57    //
58
59    return ((p->pos[Z]-Bottom()) > -kTolerance && (p->pos[Z]-Top()) < kTolerance &&
60            (p->pos.Perp()-Radius()) < kTolerance);
61}
62
63//______________________________________________________________________________
64Double_t DetectorPhotonTransporter::CylinderIntersection( const TVector3 &pos,
65       const TVector3 &dir, Bool_t zero) const {
66    //
67    // Calculate the next intersection of p with this transporter boundaries
68    //
69
70    return CylinderIntersection( pos, dir, fR, fPos[Z]+fDZup, fPos[Z]-fDZdown, zero);
71}
72
73//______________________________________________________________________________
74Double_t DetectorPhotonTransporter::CylinderIntersection( const Photon *p,
75       Bool_t zero) const {
76    //
77    // Calculate the next intersection of p with this transporter boundaries
78    //
79
80    return CylinderIntersection( p, fR, fPos[Z]+fDZup, fPos[Z]-fDZdown, zero );
81}
82
83//______________________________________________________________________________
84Double_t DetectorPhotonTransporter::CylinderIntersection( const Photon *p,
85    Double_t radius, Double_t  zup, Double_t zdown, Bool_t zero ) const {
86    //
87    // Finds, if exists, the next interaction of p over cyl.  In any case saves
88    // the distances between p->pos and the two int points in fLastCylInt
89    //
90
91    return CylinderIntersection(p->pos, p->dir, radius, zup, zdown, zero );
92}
93
94//______________________________________________________________________________
95Double_t DetectorPhotonTransporter::CylinderIntersection( const TVector3 &pos,
96    const TVector3 &dir, Double_t radius, Double_t  zup, Double_t zdown,
97    Bool_t zero ) const {
98    //
99    // Finds, if exists, the next interaction of p over cyl.  In any case saves
100    // the distances between p->pos and the two int points in fLastCylInt
101    //
102
103    Double_t dist[4] = {-kHuge,-kHuge,-kHuge,-kHuge};
104    // array of intersections, 0, 1 basis, 3, 4 lateral surface
105
106    // int with the basis
107    if ( Abs(dir[Z]) > kTolerance ) {
108        dist[0] = (zup-pos[Z])/dir[Z];
109        dist[1] = (zdown-pos[Z])/dir[Z];
110
111        //check if this points belong to the surface
112        TVector3 pos1=pos+dir*dist[0];
113        if ( pos1.Perp2()>radius*radius-kTolerance ) dist[0]=-kHuge;
114        pos1=pos+dir*dist[1];
115        if ( pos1.Perp2()>radius*radius-kTolerance ) dist[1]=-kHuge;
116    }
117
118    Double_t a, b, c;
119    a = dir[X]*dir[X]+dir[Y]*dir[Y];
120    b = 2*(dir[X]*pos[X]+dir[Y]*pos[Y]);
121    c = pos[X]*pos[X]+pos[Y]*pos[Y]-radius*radius;
122
123
124    Double_t delta = (b*b)-4*a*c;
125
126    if ( a > 0 ) {
127        if ( Abs(delta) < kTolerance ) {
128            dist[2] = -b/(2*a);
129        } else if ( delta >= kTolerance ) {
130            dist[2] = (-b-Sqrt(delta))/(2*a);
131            dist[3] = (-b+Sqrt(delta))/(2*a);
132        }
133        double newz=pos.Z()+dir.Z()*dist[2];
134        double zup_tol=zup-kTolerance;
135        double zdown_tol=zdown+kTolerance;
136        if ( (newz-zup_tol)*(newz-zdown_tol)>0 ) dist[2]=-kHuge;
137        newz=pos.Z()+dir.Z()*dist[3];
138        if ( (newz-zup_tol)*(newz-zdown_tol)>0 ) dist[3]=-kHuge;
139    }
140
141
142    Double_t threshold = zero ? -kTolerance : kTolerance;
143    Double_t dt = MaxElement(4,dist);
144    for(Int_t i(0); i<4; i++) {
145        if ( dist[i] > threshold && dist[i] < dt ) dt = dist[i];
146    }
147
148    // inbetween +/- ktolerance dt is basically 0
149    if ( zero && dt > -kTolerance && dt < kTolerance) dt = 0;
150    //MsgForm(EsafMsg::Info,"dt %3f\t,d(%8f,%8f,%8f,%8f)",dt,dist[0],dist[1],dist[2],dist[3]);
151    return dt;
152
153/***************************************************************************
154 *
155 * Old code, kept for the time being as reference
156 *
157 * *************************************************************************
158
159    // intesection point
160    TVector3 intPoint(0,0,0), uDir;
161    uDir = dir.Unit();
162
163    Double_t a, b, c, dup, ddown, dummy[2];
164
165    fLastCylInt.first=0;
166    fLastCylInt.second[0]=fLastCylInt.second[1]=0;
167
168    // distances to the upper and lower planes
169    if ( TMath::Abs(uDir[Z]) >= kTolerance ){
170        dup = (zUp-pos[Z])/uDir[Z];
171        ddown = (zDown-pos[Z])/uDir[Z];
172
173    }
174
175    if((TMath::Abs(uDir.Theta()) < kTolerance) ||
176            (TMath::Abs(uDir.Theta() - TMath::Pi()) < kTolerance)) {
177        // photon directed along Z: special case
178        // doesn't hit the cyl.
179        if(pos.Perp()-radius > kTolerance) {
180            return -1*kHuge;
181        }
182        // intersection with the bases
183        if(uDir[Z] > 0 && pos[Z] < zDown ||
184                uDir[Z] < 0 && pos[Z] > zDown && pos[Z] < zUp ){
185            // lower base
186            fLastCylInt.second[0]= ddown;
187            fLastCylInt.second[1]= dup;
188        }
189        else if(uDir[Z] < 0 && pos[Z] > zUp ||
190                uDir[Z] > 0 && pos[Z] > zDown && pos[Z] < zUp ){
191            // upper base
192            fLastCylInt.second[0]= dup;
193            fLastCylInt.second[1]= ddown;
194        } else {
195            MsgForm(EsafMsg::Warning,"CylinderIntersection(): photon lost");
196            return -2*kHuge;
197        }
198        return fLastCylInt.second[0];
199    }
200
201    // fill the coefficients of the equation
202    a=uDir[X]*uDir[X]+uDir[Y]*uDir[Y];
203    b=2*(uDir[X]*pos[X]+uDir[Y]*pos[Y]);
204    c=pos[X]*pos[X]+pos[Y]*pos[Y]-radius*radius;
205    // find the first intersection with the cylinder if it exists
206    pair<int, double* > res;
207    res=findRoots(a,b,c);
208
209    if (res.first == 0) {
210        // photon doesn't cross the cylinder
211        // return something;
212        return -3*kHuge;
213
214    } else if (res.first == 1) {
215        // one solution, tangent photon
216        intPoint=pos+uDir*res.second[0];
217        // check the sol to be inside cylinder
218        if ( intPoint[Z] > zUp || intPoint[Z] < zDown )
219            return -4*kHuge;
220        fLastCylInt.first=1;
221        fLastCylInt.second[0]=fLastCylInt.second[1]=res.second[0];
222        return fLastCylInt.second[0];
223
224    } else if ( TMath::Abs(uDir[Z]) < kTolerance ){
225        // two crossings but no intersection with the bases
226        // horizontal photon
227        intPoint=pos+uDir*res.second[0];
228        // check the sol to be inside cylinder
229        if ( intPoint[Z] > zUp || intPoint[Z] < zDown )
230            return -5*kHuge;
231
232    }
233
234    // 4 candidates left, just 2 of them on the borders of the cyl; in fact 3 if
235    // the photons crosses the intersection between the side surface and one of
236    // the bases
237    fLastCylInt.first=2;
238
239    // distances of the 4 intpoints from pos
240    // 0: upper base plane
241    // 1: lower base plane
242    // 2: 1st lateral surface
243    // 3: 2nd lateral surface
244    Double_t intDist[4];
245    intDist[0]=dup;
246    intDist[1]=ddown;
247    intDist[2]=res.second[0];
248    intDist[3]=res.second[1];
249
250    Int_t nTrueInts(0);
251    Bool_t onSide, onBase;
252    for(Int_t i(0); i<4; i++ ){
253        intPoint = pos+uDir*intDist[i];
254        onSide = TMath::Abs(intPoint.Perp()-radius) < kTolerance &&
255            (intPoint[Z] < zUp && intPoint[Z] > zDown);
256        onBase = (TMath::Abs(intPoint[Z]-zUp) < kTolerance ||
257                TMath::Abs(intPoint[Z]-zDown) < kTolerance) &&
258            intPoint.Perp() <= radius;
259        if ( onSide || onBase )
260            // if intDist[i] < kTolerance, save just 0
261            dummy[nTrueInts++]=( TMath::Abs(intDist[i]) > kTolerance ? intDist[i] : 0 );
262    }
263
264    if ( nTrueInts == 3 ) {
265        // check if two of the ints are the same point (one with the bases and
266        // one with the side)
267        if ( TMath::Abs(intDist[0]-intDist[2]) < kTolerance ||
268                TMath::Abs(intDist[0]-intDist[3]) < kTolerance ||
269                TMath::Abs(intDist[1]-intDist[2]) < kTolerance ||
270                TMath::Abs(intDist[1]-intDist[3]) < kTolerance) {
271            nTrueInts--;
272        }
273
274    }
275
276    if ( nTrueInts == 2) {
277        // last thing to do, sort the sols
278        if ( dummy[0]*dummy[1] < 0 || dummy[0]+dummy[1] < 0 ){
279            // opposite sign sols or sols < 0
280            fLastCylInt.second[0]=TMath::Max(dummy[0],dummy[1]);
281            fLastCylInt.second[1]=TMath::Min(dummy[0],dummy[1]);
282        } else {
283            // sols >= 0
284            fLastCylInt.second[0]=TMath::Min(dummy[0],dummy[1]);
285            fLastCylInt.second[1]=TMath::Max(dummy[0],dummy[1]);
286        }
287        if ( fLastCylInt.second[0] < 0 )
288            return -6*kHuge;
289        else
290            return fLastCylInt.second[0];
291
292    } else if ( nTrueInts==0 ) {
293        // non of the sols was on the cyl
294        return -7*kHuge;
295    } else {
296        // dump to screen all the intersections
297        Msg(EsafMsg::Warning)  << "pos" << pos << "   dir" << uDir << endl;
298        for(Int_t i(0); i<4; i++ ){
299            MsgForm(EsafMsg::Warning, "Intersection %d:",i);
300            intPoint = pos+uDir*intDist[i];
301            MsgForm(EsafMsg::Warning,"intPoint.PerP() = %f intPoint[Z] = %f",
302                    intPoint.Perp(),intPoint[Z]);
303            MsgForm(EsafMsg::Warning,"Radius = %f Zup = %f Zdown = %f",
304                    radius, zUp, zDown);
305
306            onSide = TMath::Abs(intPoint.Perp()-radius) < kTolerance &&
307                (intPoint[Z] < zUp && intPoint[Z] > zDown);
308
309            onBase = (TMath::Abs(intPoint[Z]-zUp) < kTolerance ||
310                    TMath::Abs(intPoint[Z]-zDown) < kTolerance) &&
311                intPoint.Perp()-radius <= kTolerance;
312            MsgForm(EsafMsg::Warning," side = %f base = %f", onSide, onBase);
313        }
314        Msg(EsafMsg::Warning) << " nTrueInts = " << nTrueInts << MsgDispatch;
315        MsgForm(EsafMsg::Panic,"CylinderIntersection(): something gone wrong!");
316        return -8*kHuge;
317    }
318*******************************************************************************/
319}
Note: See TracBrowser for help on using the repository browser.