Corrected drawing routine for RecHit3D (J. Barbosa)
[u/mrichter/AliRoot.git] / RICH / AliRICHEllipse.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17   $Log$
18   Revision 1.4  2001/05/10 12:34:43  jbarbosa
19   Changed drwaing routines.
20
21   Revision 1.3  2000/11/01 15:37:44  jbarbosa
22   Removed verbose output.
23
24   Revision 1.2  2000/06/30 16:31:51  dibari
25   New drawing routine from Nico and Daniela.
26
27   Revision 1.1  2000/06/12 15:21:57  jbarbosa
28   Cleaned up version.
29
30 */
31
32 #include "TMath.h"
33 #include "AliRICHEllipse.h"
34 #include "AliRICH.h"
35 #include "AliRun.h"
36 #include "AliRICHPatRec.h"
37 #include "AliRICHGeometry.h"
38
39 #include <TRandom.h>
40
41 ClassImp(AliRICHEllipse)
42
43 //________________________________________________________________________________
44 AliRICHEllipse::AliRICHEllipse()
45
46
47 //  Default Constructor for a RICH ellipse
48
49     fCx = 0;
50     fCy = 0;
51     fOmega = 0;
52     fTheta = 0;
53     fPhi = 0;
54     fh= 0;
55 }
56
57 //________________________________________________________________________________
58 AliRICHEllipse::~AliRICHEllipse()
59
60
61 // Destructor
62
63     fCx = 0;
64     fCy = 0;
65     fOmega = 0;
66     fTheta = 0;
67     fPhi = 0;
68     fh= 0;
69 }
70
71
72 //________________________________________________________________________________
73 AliRICHEllipse::AliRICHEllipse(Float_t cx, Float_t cy, Float_t omega, Float_t theta, Float_t phi, Float_t emiss)
74
75
76 //  Constructor for a RICH ellipse
77
78     fCx = cx;
79     fCy = cy;
80     fOmega = omega;
81     fTheta = theta;
82     fPhi = phi;
83     fEmissPoint = emiss;
84     fh=9.25;
85 }
86
87 //________________________________________________________________________________
88 void AliRICHEllipse::CreatePoints(Int_t chamber)
89 {
90
91 // Create points along the ellipse equation
92   
93   Float_t x, y, rotx, roty, h, cx, cy, phi, omega, theta, omega1, theta1, phiinc;
94   Float_t a,b,c,r,e, offset;
95   
96   Float_t kPi=TMath::Pi();
97
98   AliRICH *pRICH  = (AliRICH*)gAlice->GetModule("RICH");
99   AliRICHChamber*       iChamber;
100   AliRICHGeometry*  geometry;
101   
102   iChamber = &(pRICH->Chamber(chamber));
103   geometry=iChamber->GetGeometryModel();
104   
105   //h = 2.3 * geometry->GetRadiatorToPads();
106   h = geometry->GetRadiatorToPads();
107   //printf("h: %f",h);
108
109   cx = fCx;
110   cy = fCy;
111   theta = fTheta;
112   omega = fOmega;
113   phiinc = fPhi+kPi/2;
114   
115   printf("Omega: %f, Theta: %f, Phi: %f\n", omega, theta, phiinc); 
116
117
118   for(Float_t i=0;i<1000;i++)
119     {
120       phi=((Float_t)(i)/1000)*2*kPi;
121       //printf("Phi: %f\n",phi);
122
123       //theta1=SnellAngle(theta1);
124       
125       //if(phi<=TMath::Pi())
126       omega1=SnellAngle(omega);
127       theta1=SnellAngle(theta);
128       //omega1=SnellAngle(omega+cos(phi)*theta);
129       //if(phi>TMath::Pi())
130         //omega1=SnellAngle(omega+(1-2*(2*TMath::Pi()-phi)/(TMath::Pi()))*theta);
131
132    
133       //Omega1->Fill(omega1,(float) 1);
134
135       a = h*(tan(omega1+theta1)+tan(omega1-theta1))/2;
136       b = h*tan(omega1);
137       e = sqrt(1 - (b*b)/(a*a));
138       c = a*e;
139       r = (a*(1-e*e))/(1+e*cos(e));
140       offset = h*(tan(omega1+theta1)-tan(omega1-theta1))/2;
141         
142       x = b* sin(phi);
143       y = a* cos(phi) + offset;
144                 
145       rotx = x*cos(phiinc)-y*sin(phiinc);
146       roty = x*sin(phiinc)+y*cos(phiinc);
147    
148       //x = h * 1/(tan(omega1)) * sin(phi+phiinc);
149       //y = x * 1/(tan(phi+phiinc));
150
151       
152
153       //Rings->Fill(x,y,(float) 1);
154
155       rotx += cx;
156       roty += cy;
157
158       //printf("x:%f, y: %f\n",x,y);
159
160       Float_t vectorLoc[3]={rotx,6.276,roty};
161       Float_t  vectorGlob[3];
162       iChamber->LocaltoGlobal(vectorLoc,vectorGlob);
163       SetPoint((Int_t) i,vectorGlob[0],vectorGlob[1],vectorGlob[2]);
164       //printf("Coordinates: %f %f %f\n",vectorGlob[0],vectorGlob[1],vectorGlob[2]);
165       
166     }
167
168 }
169
170 void AliRICHEllipse::CerenkovRingDrawing(Int_t chamber,Int_t track)
171 {
172
173 //to draw Cherenkov ring by known Cherenkov angle
174
175   Int_t nmaxdegrees;
176   Int_t Nphpad;
177   Float_t phpad;
178   Float_t nfreonave, nquartzave;
179   Float_t aveEnerg;
180   Float_t energy[2];
181   Float_t e1, e2, f1, f2;
182   Float_t pointsOnCathode[3];
183
184   //printf("Drawing ring in chamber:%d\n",chamber);
185
186
187   AliRICH *pRICH  = (AliRICH*)gAlice->GetModule("RICH");
188   AliRICHChamber*       iChamber;
189   
190   iChamber = &(pRICH->Chamber(chamber));
191
192   AliRICHPatRec *PatRec = new AliRICHPatRec;
193   PatRec->TrackParam(track,chamber,fTheta,fOmega);
194
195   //printf("Just created PateRec\n");
196
197 //parameters to calculate freon window refractive index vs. energy
198
199     Float_t a = 1.177;
200     Float_t b = 0.0172;
201     
202 //parameters to calculate quartz window refractive index vs. energy
203 /*
204    Energ[0]  = 5.6;
205    Energ[1]  = 7.7;
206 */      
207     energy[0]  = 5.0;
208     energy[1]  = 8.0;
209     e1  = 10.666;
210     e2  = 18.125;
211     f1  = 46.411;
212     f2  = 228.71;
213   
214
215     /*Float_t nquartz = 1.585;
216       Float_t ngas    = 1.;
217       Float_t nfreon  = 1.295;
218       Float_t value;
219     */
220
221
222
223    nmaxdegrees = 360;
224    
225    for (Nphpad=0; Nphpad<nmaxdegrees;Nphpad++) { 
226
227        phpad = (360./(Float_t)nmaxdegrees)*(Float_t)Nphpad;
228       
229        aveEnerg =  (energy[0]+energy[1])/2.;
230        //aveEnerg = 6.5;
231        
232        
233        nfreonave  = a+b*aveEnerg;
234        nquartzave = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(aveEnerg,2)))+
235                          (f2/(TMath::Power(e2,2)-TMath::Power(aveEnerg,2))));
236
237        //nfreonave = 1.295;
238        //nquartzave = 1.585;
239        
240        ///printf("Calling DistancefromMip %f %f \n",fEmissPoint,fOmega);
241        
242        //Float_t dummy = 
243          PatRec->DistanceFromMip(nfreonave, nquartzave,fEmissPoint,fOmega, phpad, pointsOnCathode,fTheta,fPhi);
244
245        //Float_t points[3];
246
247        //points = pointsOnCathode;
248
249
250        //printf(" values %f %f %f\n",points[0],points[1],points[2]);
251        
252        Float_t vectorLoc[3]={pointsOnCathode[0],pointsOnCathode[2],pointsOnCathode[1]};
253        Float_t  vectorGlob[3];
254        iChamber->LocaltoGlobal(vectorLoc,vectorGlob);
255        SetPoint(Nphpad,vectorGlob[0],vectorGlob[1],vectorGlob[2]);
256       //fCoordEllipse[0][Nphpad] = pointsOnCathode[0];
257       //fCoordEllipse[1][Nphpad] = pointsOnCathode[1];
258        
259        //printf(" values %f %f \n",pointsOnCathode[0],pointsOnCathode[1]);
260        
261    }
262
263 }
264
265
266
267 Float_t AliRICHEllipse:: SnellAngle(Float_t iangle)
268
269
270 // Compute the Snell angle
271
272   Float_t nfreon  = 1.295;
273   Float_t nquartz = 1.585;
274   Float_t ngas = 1;
275
276   Float_t sinrangle;
277   Float_t rangle;
278   Float_t a1, a2;
279
280   a1=nfreon/nquartz;
281   a2=nquartz/ngas;
282
283   sinrangle = a1*a2*sin(iangle);
284
285   if(sinrangle>1.) {
286      rangle = 999.;
287      return rangle;
288   }
289   
290   rangle = asin(sinrangle);  
291   //printf("iangle %f, a1*a2, %f, sinranlge, %f, rangle, %f\n", iangle, a1*a2, sinrangle, rangle);
292   return rangle;
293
294 }
295
296
297
298
299
300
301
302