Cleaned up version.
[u/mrichter/AliRoot.git] / RICH / AliRICHPatRec.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.1  2000/06/09 14:53:01  jbarbosa
19   Bari's pattern recognition algorithm
20
21 */
22
23 #include "AliRICHHit.h"
24 #include "AliRICHCerenkov.h"
25 #include "AliRICHPadHit.h"
26 #include "AliRICHDigit.h"
27 #include "AliRICHRawCluster.h"
28 #include "AliRICHRecHit.h"
29 #include "AliRun.h"
30 #include "AliDetector.h"
31 #include "AliRICH.h"
32 #include "AliRICHPoints.h"
33 #include "AliRICHSegmentation.h"
34 #include "AliRICHPatRec.h"
35 #include "AliRICH.h"
36 #include "AliRICHConst.h"
37 #include "AliRICHPoints.h"
38 #include "AliConst.h"
39
40 #include <TParticle.h>
41 #include <TMath.h>
42 #include <TRandom.h>
43 #include <TCanvas.h>
44 #include <TH2.h>
45
46
47 ClassImp(AliRICHPatRec)
48 //___________________________________________
49 AliRICHPatRec::AliRICHPatRec() : TObject()
50 {
51   // Default constructor
52   
53     //fChambers = 0;
54 }
55 //___________________________________________
56 AliRICHPatRec::AliRICHPatRec(const char *name, const char *title)
57     : TObject()
58 {
59     //Constructor for Bari's pattern recogniton method object
60 }
61
62 void AliRICHPatRec::PatRec()
63 {
64
65 // Pattern recognition algorithm
66
67   AliRICHChamber*       iChamber;
68   AliRICHSegmentation*  segmentation;
69         
70   Int_t ntracks, ndigits[kNCH];
71   Int_t itr, ich, i;
72   Int_t goodPhotons;
73   Int_t x,y,q;
74   Float_t rx,ry;
75   Int_t nent,status;
76
77   Float_t gamma,massCer,betaCer;
78
79   Float_t rechit[5];
80
81   printf("PatRec started\n");
82
83   TCanvas *c1    = new TCanvas("c1","Alice RICH pad hits",50,10,700,700);
84
85   TH2F *ring     = new TH2F("ring","map",90,-30.,30.,90,-30.,30.);
86   TH2F *ringband = new TH2F("ringband","map",90,-65.,65.,90,-65.,65.);
87   TH1F *cerangle = new TH1F("cerangle","phot",300,0.45,0.75);
88   TH1F *ceranglew= new TH1F("ceranglew","phot",300,0.45,0.75);
89   TH1F *hough    = new TH1F("hough","hough",75,0.45,0.75);
90   TH1F *mass     = new TH1F("mass","mass",100,50.,1050.); 
91  
92   AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
93   TTree *treeH = gAlice->TreeH();
94
95   ntracks =(Int_t) treeH->GetEntries();
96   //  ntracks = 1;
97   for (itr=0; itr<ntracks; itr++) {
98  
99     status = TrackParam(itr,ich);    
100     if(status==1) continue;
101     printf(" theta %f phi %f track \n",fTrackTheta,fTrackPhi);
102     //    ring->Fill(fTrackLoc[0],fTrackLoc[1],100.);
103
104     iChamber = &(pRICH->Chamber(ich));
105     segmentation=iChamber->GetSegmentationModel();
106
107     nent=(Int_t)gAlice->TreeD()->GetEntries();
108     gAlice->TreeD()->GetEvent(nent-1);
109     TClonesArray *pDigitss = pRICH->DigitsAddress(ich);
110     ndigits[ich] = pDigitss->GetEntriesFast();
111     printf("ndigits %d in chamber %d\n",ndigits[ich],ich);
112     AliRICHDigit *padI = 0;
113
114     goodPhotons = 0;
115
116     for (Int_t dig=0;dig<ndigits[ich];dig++) {
117       padI=(AliRICHDigit*) pDigitss->UncheckedAt(dig);
118       x=padI->fPadX;
119       y=padI->fPadY;
120       q=padI->fSignal;
121       segmentation->GetPadCxy(x,y,rx,ry);      
122
123       fXpad = rx-fXshift;
124       fYpad = ry-fYshift;
125       fQpad = q;
126     
127       ringband->Fill(x,y,1.);
128       fCerenkovAnglePad = PhotonCerenkovAngle();
129       if(fCerenkovAnglePad==-999) continue;
130
131       if(!PhotonInBand()) continue;
132
133       ring->Fill(fXpad,fYpad,1.);
134       cerangle->Fill(fCerenkovAnglePad,1.);
135
136       goodPhotons++;
137       fEtaPhotons[goodPhotons] = fCerenkovAnglePad;
138     }
139     fNumEtaPhotons = goodPhotons;
140
141     BackgroundEstimation();
142
143     for(i=0;i<goodPhotons;i++) {
144       ceranglew->Fill(fEtaPhotons[i],fWeightPhotons[i]);
145       //      printf(" Eta %f weight %f \n",fEtaPhotons[i],fWeightPhotons[i]);
146     }
147
148     HoughResponse();
149
150     rechit[0] = 0;
151     rechit[1]   = 0;
152     rechit[2] = fThetaCerenkov;
153     rechit[3] = 0;
154     rechit[4] = 0;
155
156
157     hough->Fill(fThetaCerenkov,1.);
158     
159     pRICH->AddRecHit(ich,rechit);
160     
161     betaCer = BetaCerenkov(1.29,fThetaCerenkov);
162     gamma  = 1./sqrt(1.-pow(betaCer,2));
163     massCer = fTrackMom/(betaCer*gamma);
164     //    printf(" mass %f \n",massCer);
165     mass->Fill(massCer*1000,1.);
166   }    
167
168   gAlice->TreeR()->Fill();
169   TClonesArray *fRec;
170   for (i=0;i<kNCH;i++) {
171     fRec=pRICH->RecHitsAddress(i);
172     int ndig=fRec->GetEntriesFast();
173     printf ("Chamber %d, rings %d\n",i,ndig);
174   }
175   pRICH->ResetRecHits();
176   
177
178   c1->Divide(2,2);
179   c1->cd(1);
180   ring->Draw("box");
181   c1->cd(2);
182   ringband->Draw("box");
183   c1->cd(3);
184   cerangle->Draw();
185   c1->cd(4);
186   hough->Draw();
187
188 }     
189
190
191 Int_t AliRICHPatRec::TrackParam(Int_t itr, Int_t &ich)
192 {
193   // Get Local coordinates of track impact  
194
195   AliRICHChamber*       iChamber;
196   AliRICHSegmentation*  segmentation;
197
198   Float_t trackglob[3];
199   Float_t trackloc[3];
200   Float_t thetatr;  
201   Float_t phitr;
202   Float_t iloss;
203   Float_t part;
204   Float_t pX, pY, pZ;
205
206   printf("Calling TrackParam\n");
207
208     gAlice->ResetHits();
209     TTree *treeH = gAlice->TreeH();
210     treeH->GetEvent(itr);
211  
212     AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
213     AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
214     if(mHit==0) return 1;
215     ich = mHit->fChamber-1;
216     trackglob[0] = mHit->fX;
217     trackglob[1] = mHit->fY;
218     trackglob[2] = mHit->fZ;
219     pX = mHit->fMomX;
220     pY = mHit->fMomY;
221     pZ = mHit->fMomZ;
222     fTrackMom = sqrt(pow(pX,2)+pow(pY,2)+pow(pZ,2));
223     thetatr = (180 - mHit->fTheta)*(Float_t)kDegrad;
224     phitr = mHit->fPhi*(Float_t)kDegrad;
225     iloss = mHit->fLoss;
226     part  = mHit->fParticle;
227
228     iChamber = &(pRICH->Chamber(ich));
229     iChamber->GlobaltoLocal(trackglob,trackloc);
230
231     segmentation=iChamber->GetSegmentationModel();
232
233     // retrieve geometrical params
234
235     AliRICHGeometry* fGeometry=iChamber->GetGeometryModel();   
236     
237     fRw   = fGeometry->GetFreonThickness();
238     fQw   = fGeometry->GetQuartzThickness();
239     fTgap = fGeometry->GetGapThickness() 
240             + fGeometry->GetProximityGapThickness();  
241
242     Float_t apar = (fRw + fQw + fTgap)*tan(thetatr);  
243     fTrackLoc[0] = apar*cos(phitr);
244     fTrackLoc[1] = apar*sin(phitr);
245     fTrackLoc[2] = fRw + fQw + fTgap;
246     fTrackTheta  = thetatr;
247     fTrackPhi    = phitr;
248    
249     fXshift = trackloc[0] - fTrackLoc[0];    
250     fYshift = trackloc[2] - fTrackLoc[1];
251      
252     return 0;
253 }
254
255 Float_t AliRICHPatRec::EstimationAtLimits(Float_t lim, Float_t radius,
256                                                  Float_t phiphot)
257 {
258
259 // Estimation of emission point
260
261   Float_t nquartz = 1.585;
262   Float_t ngas    = 1.;
263   Float_t nfreon  = 1.295;
264   Float_t value;
265
266   //  printf("Calling EstimationLimits\n");
267
268   Float_t apar = (fRw -fEmissPoint + fQw + fTgap)*tan(fTrackTheta);
269   Float_t b1 = (fRw-fEmissPoint)*tan(lim);
270   Float_t b2 = fQw / sqrt(pow(nquartz,2)-pow(nfreon*sin(lim),2));
271   Float_t b3 = fTgap / sqrt(pow(ngas,2)-pow(nfreon*sin(lim),2));
272   Float_t bpar = b1 + nfreon*sin(lim)*(b2+b3);
273   value = pow(radius,2)
274     -pow((apar*cos(fTrackPhi)-bpar*cos(phiphot)),2)
275     -pow((apar*sin(fTrackPhi)-bpar*sin(phiphot)),2);
276   return value;
277 }
278
279
280 Float_t AliRICHPatRec::PhotonCerenkovAngle()
281 {
282   // Cherenkov pad angle reconstruction
283      
284   Float_t radius;
285   Float_t cherMin = 0;
286   Float_t cherMax = 0.8;
287   Float_t phiphot;
288   Float_t eps = 0.0001;
289   Int_t   niterEmiss = 0;
290   Int_t   niterEmissMax = 0;
291   Float_t x1,x2,x3=0,p1,p2,p3;
292   Float_t argY,argX;
293   Int_t niterFun;
294
295   //  printf("Calling PhotonCerenkovAngle\n");
296
297   radius = sqrt(pow(fTrackLoc[0]-fXpad,2)+pow(fTrackLoc[1]-fYpad,2));
298   fEmissPoint = fRw/2.;  //Start value of EmissionPoint
299   
300   while(niterEmiss<=niterEmissMax) {
301  
302     niterFun = 0;
303     argY = fYpad - fEmissPoint*tan(fTrackTheta)*sin(fTrackPhi);
304     argX = fXpad - fEmissPoint*tan(fTrackTheta)*cos(fTrackPhi);
305     phiphot = atan2(argY,argX); 
306     p1 = EstimationAtLimits(cherMin,radius,phiphot);
307     p2 = EstimationAtLimits(cherMax,radius,phiphot);
308     if(p1*p2>0)
309       {
310         //        printf("PhotonCerenkovAngle failed\n");
311         return -999;
312       }
313
314     //start to find the Cherenkov pad angle 
315     x1 = cherMin;
316     x2 = cherMax;
317     x3 = (x1+x2)/2.;
318     p3 = EstimationAtLimits(x3,radius,phiphot);
319     while(TMath::Abs(p3)>eps){
320       if(p1*p3<0) x2 = x3;
321       if(p1*p3>0) {
322         x1 = x3;
323         p1 = EstimationAtLimits(x1,radius,phiphot);
324       }
325       x3 = (x1+x2)/2.;
326       p3 = EstimationAtLimits(x3,radius,phiphot);
327       niterFun++;
328
329       if(niterFun>=1000) {
330         //      printf(" max iterations in PhotonCerenkovAngle\n");
331         return x3;
332       }
333     }
334     //    printf("niterFun %i \n",niterFun);
335     niterEmiss++;
336     if (niterEmiss != niterEmissMax+1) EmissionPoint();
337   } 
338   /* 
339    printf(" phiphot %f fXpad %f fYpad %f fEmiss %f \n",
340                          phiphot,fXpad,fYpad,fEmissPoint);
341   */
342
343   return x3;
344
345 }
346
347
348 void AliRICHPatRec::EmissionPoint()
349
350
351 // Find emission point
352
353   Float_t absorbtionLength=7.83*fRw; //absorption length in the freon (cm)
354   // 7.83 = -1/ln(T0) where 
355   // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
356   Float_t photonLength, photonLengthMin, photonLengthMax;
357
358   photonLength=exp(-fRw/(absorbtionLength*cos(fCerenkovAnglePad)));
359   photonLengthMin=fRw*photonLength/(1.-photonLength);
360   photonLengthMax=absorbtionLength*cos(fCerenkovAnglePad);
361   fEmissPoint = fRw + photonLengthMin - photonLengthMax;
362
363 }
364
365 void AliRICHPatRec::PhotonSelection(Int_t track, Int_t &nphot, Float_t &thetamean)
366 {
367
368 // not implemented yet
369
370   printf("Calling PhotonSelection\n");
371 }
372
373 void AliRICHPatRec::BackgroundEstimation()
374 {
375
376 // estimate background noise
377
378   Float_t stepEta   = 0.001;  
379   Float_t etaMinBkg = 0.72;
380   Float_t etaMaxBkg = 0.75;
381   Float_t etaMin    = 0.;
382   Float_t etaMax    = 0.75;
383   Float_t ngas      = 1.;
384   Float_t nfreon    = 1.295;
385
386   Float_t etaStepMin,etaStepMax,etaStepAvg;
387   Int_t i,ip,nstep;
388   Int_t numPhotBkg, numPhotonStep;
389   Float_t funBkg,areaBkg,normBkg;
390   Float_t densityBkg,storeBkg,numStore;
391   Float_t thetaSig;
392   
393   numPhotBkg = 0;
394   areaBkg = 0.;
395
396   nstep = (int)((etaMaxBkg-etaMinBkg)/stepEta);
397
398   for (i=0;i<fNumEtaPhotons;i++) {
399
400     if(fEtaPhotons[i]>etaMinBkg && fEtaPhotons[i]<etaMaxBkg) {
401       numPhotBkg++;
402     }    
403   }
404   if (numPhotBkg == 0) {
405      for (i=0;i<fNumEtaPhotons;i++) {
406         fWeightPhotons[i] = 1.;
407      }
408     return;
409   }
410
411   //  printf(" numPhotBkg %i ",numPhotBkg);
412
413   for (i=0;i<nstep;i++) {
414     etaStepMin = etaMinBkg + (Float_t)(i)*stepEta;
415     etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;    
416     etaStepAvg = 0.5*(etaStepMax + etaStepMin);
417     /*
418     funBkg = tan(etaStepAvg)*pow((1.+pow(tan(etaStepAvg),2)),
419                                   5.52)-7.803 + 22.02*tan(etaStepAvg);
420     */
421     thetaSig = asin(nfreon/ngas*sin(etaStepAvg));
422     funBkg = tan(thetaSig)*(1.+pow(tan(thetaSig),2))*nfreon
423        /ngas*cos(etaStepAvg)/cos(thetaSig);
424     areaBkg += stepEta*funBkg;
425   }
426
427   densityBkg = 0.95*(Float_t)(numPhotBkg)/areaBkg;
428   //  printf(" densityBkg %f \n",densityBkg);
429   
430   nstep = (int)((etaMax-etaMin)/stepEta); 
431   storeBkg = 0.;
432   numStore = 0;
433   for (i=0;i<nstep;i++) {
434     etaStepMin = etaMinBkg + (Float_t)(i)*stepEta;
435     etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;    
436     etaStepAvg = 0.5*(etaStepMax + etaStepMin);
437     /*
438     funBkg = tan(etaStepAvg)*pow((1.+pow(tan(etaStepAvg),2)),
439                                   5.52)-7.803 + 22.02*tan(etaStepAvg);
440     */
441
442     thetaSig = asin(nfreon/ngas*sin(etaStepAvg));
443     funBkg = tan(thetaSig)*(1.+pow(tan(thetaSig),2))*nfreon
444        /ngas*cos(etaStepAvg)/cos(thetaSig);
445
446     areaBkg = stepEta*funBkg;
447     normBkg = densityBkg*areaBkg; 
448     numPhotonStep = 0;
449     for (ip=0;ip<fNumEtaPhotons;ip++) {
450       if(fEtaPhotons[ip]>etaStepMin && fEtaPhotons[ip]<etaStepMax) {
451         numPhotonStep++;
452       }
453     }
454     if (numPhotonStep == 0) {
455       storeBkg += normBkg;
456       numStore++;
457       if (numStore>50) {
458         numStore = 0;
459         storeBkg = 0.;
460       }
461     }
462     if (numPhotonStep == 0) continue;
463     for (ip=0;ip<fNumEtaPhotons;ip++) {
464       if(fEtaPhotons[ip]>etaStepMin && fEtaPhotons[ip]<etaStepMax) {
465         normBkg +=storeBkg;
466         storeBkg = 0;
467         numStore = 0;
468         fWeightPhotons[ip] = 1. - normBkg/(Float_t)(numPhotonStep);
469         /*
470         printf(" normBkg %f numPhotonStep %i fW %f \n",
471                normBkg, numPhotonStep, fWeightPhotons[ip]);
472         */
473         if(fWeightPhotons[ip]<0) fWeightPhotons[ip] = 0.;
474       }
475     }
476   }
477 }
478
479
480 void AliRICHPatRec::FlagPhotons(Int_t track, Float_t theta)
481 {
482
483 // not implemented yet
484
485   printf("Calling FlagPhotons\n");
486 }
487
488
489 //////////////////////////////////////////
490
491
492
493
494
495 Int_t AliRICHPatRec::PhotonInBand()
496
497   //0=label for parameters giving internal band ellipse
498   //1=label for parameters giving external band ellipse  
499
500   Float_t imp[2], mass[2], energy[2], beta[2]; 
501   Float_t emissPointLength[2];
502   Float_t e1, e2, f1, f2;
503   Float_t nfreon[2], nquartz[2]; 
504   Int_t times;
505
506
507   Float_t phpad, thetacer[2]; 
508   Float_t bandradius[2], padradius;
509
510   imp[0] = 5.0; //threshold momentum for the proton Cherenkov emission
511   imp[1] = 1.2;
512  
513   mass[0] = 0.938; //proton mass 
514   mass[1] = 0.139; //pion mass
515
516   emissPointLength[0] = fRw-0.0001; //at the beginning of the radiator
517   emissPointLength[1] = 0.;//at the end of radiator
518   
519   //parameters to calculate freon window refractive index vs. energy
520   Float_t a = 1.177;
521   Float_t b = 0.0172;
522  
523   //parameters to calculate quartz window refractive index vs. energy
524   /*
525   Energ[0]  = 5.6;
526   Energ[1]  = 7.7;
527   */
528   energy[0]  = 5.0;
529   energy[1]  = 8.0;
530   e1 = 10.666;
531   e2  = 18.125;
532   f1  = 46.411;
533   f2  = 228.71;
534
535
536   phpad = PhiPad();  
537
538   for (times=0; times<=1; times++) {
539   
540     nfreon[times]   = a+b*energy[times];
541
542     nquartz[times] = sqrt(1+(f1/(pow(e1,2)-pow(energy[times],2)))+
543                           (f2/(pow(e2,2)-pow(energy[times],2))));
544
545     beta[times]  = imp[times]/sqrt(pow(imp[times],2)+pow(mass[times],2));
546    
547     thetacer[times] =  CherenkovAngle( nfreon[times], beta[times]);
548
549     bandradius[times] = DistanceFromMip( nfreon[times], nquartz[times],
550                                         emissPointLength[times], 
551                                         thetacer[times], phpad);
552   }
553
554   bandradius[0] -= 1.6;
555   bandradius[1] += 1.6;
556   padradius = sqrt(pow(fXpad,2)+pow(fYpad,2));
557   //  printf(" rmin %f r %f rmax %f \n",bandradius[0],padradius,bandradius[1]);
558
559   if(padradius>=bandradius[0] && padradius<=bandradius[1]) return 1;
560   return 0; 
561 }
562
563 Float_t AliRICHPatRec::DistanceFromMip(Float_t nfreon, Float_t nquartz, 
564                        Float_t emissPointLength, Float_t thetacer, 
565                        Float_t phpad)
566
567
568 // Find the distance to MIP impact
569
570   Float_t distanceValue;
571
572   TVector3 radExitPhot(1,1,1);//photon impact at the radiator exit with respect
573   //to local reference sistem with the origin in the MIP entrance 
574    
575   TVector3 vectEmissPointLength(1,1,1);
576   Float_t magEmissPointLenght;
577
578   TVector3 radExitPhot2(1,1,1);//photon impact at the radiator exit with respect
579   Float_t magRadExitPhot2;
580   //to a reference sistem with origin in the photon emission point and  
581   //axes parallel to the MIP reference sistem
582
583   TVector3 quarExitPhot(1,1,1);//photon impact at the quartz exit with respect
584   Float_t magQuarExitPhot;
585   // 
586   TVector3 gapExitPhot(1,1,1) ;
587   Float_t magGapExitPhot;
588   //
589   TVector3 fPhotocatExitPhot(1,1,1);
590   Double_t theta2;
591   Double_t thetarad , phirad ;
592   Double_t thetaquar, phiquar;
593   Double_t thetagap , phigap ;
594
595   Float_t ngas    = 1.;
596
597   magEmissPointLenght =  emissPointLength/cos(fTrackTheta);
598
599   vectEmissPointLength.SetMag(magEmissPointLenght);
600   vectEmissPointLength.SetTheta(fTrackTheta);
601   vectEmissPointLength.SetPhi(fTrackPhi);
602
603
604   radExitPhot2.SetTheta(thetacer);  
605   radExitPhot2.SetPhi(phpad); 
606
607
608   TRotation r1;
609   TRotation r2;
610   TRotation r;
611
612   r1. RotateY(fTrackTheta);
613   r2. RotateZ(fTrackPhi);
614   
615
616
617   r = r2 * r1;//rotation about the y axis by MIP theta incidence angle
618   //following by a rotation about the z axis by MIP phi incidence angle; 
619
620
621   radExitPhot2    = r * radExitPhot2;
622   theta2          = radExitPhot2.Theta();
623   magRadExitPhot2 = (fRw -  vectEmissPointLength(2))/cos(theta2);
624   radExitPhot2.SetMag(magRadExitPhot2);
625
626
627   radExitPhot = vectEmissPointLength + radExitPhot2;
628   thetarad    = radExitPhot.Theta();
629
630   phirad  =  radExitPhot.Phi(); //check on the original file //
631
632   thetaquar   = SnellAngle( nfreon, nquartz, theta2); 
633   phiquar     = radExitPhot2.Phi(); 
634   if(thetaquar == 999.) return thetaquar;
635   magQuarExitPhot    = fQw/cos(thetaquar);
636   quarExitPhot.SetMag( magQuarExitPhot);
637   quarExitPhot.SetTheta(thetaquar);
638   quarExitPhot.SetPhi(phiquar);
639
640   thetagap = SnellAngle( nquartz, ngas, thetaquar); 
641   phigap   = phiquar; 
642   if(thetagap == 999.) return thetagap;
643   magGapExitPhot    = fTgap/cos(thetagap);
644   gapExitPhot.SetMag( magGapExitPhot);
645   gapExitPhot.SetTheta(thetagap);
646   gapExitPhot.SetPhi(phigap);
647
648   fPhotocatExitPhot =  radExitPhot + quarExitPhot + gapExitPhot; 
649
650   distanceValue = sqrt(pow(fPhotocatExitPhot(0),2)
651                            +pow(fPhotocatExitPhot(1),2)); 
652   return  distanceValue ;
653 }
654
655 Float_t AliRICHPatRec::PhiPad()
656 {
657
658 // ??
659
660   Float_t zpad;
661   Float_t thetapad, phipad;
662   Float_t thetarot, phirot;
663
664   zpad = fRw + fQw + fTgap;
665
666   TVector3 photonPad(fXpad, fYpad, zpad);
667   thetapad = photonPad.Theta();
668   phipad = photonPad.Phi();
669
670   TRotation r1;
671   TRotation r2;
672   TRotation r;
673
674   thetarot = - fTrackTheta;
675   phirot   = - fTrackPhi;
676   r1. RotateZ(phirot);
677   r2. RotateY(thetarot);
678
679   r = r2 * r1;//rotation about the z axis by MIP -phi incidence angle
680   //following by a rotation about the y axis by MIP -theta incidence angle; 
681
682   photonPad  = r * photonPad;
683
684   phipad = photonPad.Phi(); 
685
686   return phipad;
687 }
688
689 Float_t AliRICHPatRec:: SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
690
691
692 // Compute the Snell angle
693
694   Float_t sinrefractangle;
695   Float_t refractangle;
696
697   sinrefractangle = (n1/n2)*sin(theta1);
698
699   if(sinrefractangle>1.) {
700      refractangle = 999.;
701      return refractangle;
702   }
703   
704   refractangle = asin(sinrefractangle);  
705   return refractangle;
706 }
707
708 Float_t AliRICHPatRec::CherenkovAngle(Float_t n, Float_t beta)
709
710
711 // Compute the cerenkov angle
712
713   Float_t thetacer;  
714       
715   if((n*beta)<1.) {
716     thetacer = 999.;
717     return thetacer;
718   }
719
720   thetacer = acos (1./(n*beta));
721   return thetacer;
722 }
723
724 Float_t AliRICHPatRec::BetaCerenkov(Float_t n, Float_t theta)
725
726
727 // Find beta
728
729   Float_t beta;  
730       
731   beta = 1./(n*cos(theta));
732   return beta;
733 }
734
735
736
737
738 void AliRICHPatRec::HoughResponse()
739
740 {       
741
742 // Implement Hough response pat. rec. method
743
744   int           bin=0;
745   int           bin1=0;
746   int           bin2=0;
747   int           i, j, k, nCorrBand;
748   int           etaBin = 750;
749   float         hcs[750];
750   float         angle, thetaCerMean;
751
752   float         etaPeak[30];
753   float         etaMin = 0.00;
754   float         etaMax = 0.75;
755   float         stepEta = 0.001;
756   float         windowEta = 0.040;
757
758   int           nBin;
759
760   float etaPeakPos  = -1;
761   Int_t   etaPeakCount = -1;
762   
763   thetaCerMean   = 0.;
764   fThetaCerenkov = 0.;    
765     
766   nBin = (int)(0.5+etaMax/(stepEta));
767   nCorrBand = (int)(0.5+ windowEta/(2 * stepEta)); 
768   memset ((void *)hcs, 0, etaBin*sizeof(int));
769
770   for (k=0; k< fNumEtaPhotons; k++) {
771
772     angle = fEtaPhotons[k];
773
774     if (angle>=etaMin && angle<= etaMax) {
775       bin = (int)(0.5+angle/(stepEta));
776       bin1= bin-nCorrBand;
777       bin2= bin+nCorrBand;
778       if (bin1<0)    bin1=0;
779       if (bin2>nBin) bin2=nBin;
780       
781       for (j=bin1; j<bin2; j++) {
782         hcs[j] += fWeightPhotons[k]; 
783       }
784
785       thetaCerMean += angle;
786     }
787   }
788  
789  thetaCerMean /= fNumEtaPhotons; 
790  
791   HoughFiltering(hcs);
792
793   for (bin=0; bin <nBin; bin++) {
794     angle = (bin+0.5) * (stepEta);
795     if (hcs[bin] && hcs[bin] > etaPeakPos) {
796       etaPeakCount = 0;
797       etaPeakPos = hcs[bin];
798       etaPeak[0]=angle;
799     }
800     else { 
801       if (hcs[bin] == etaPeakPos) {
802         etaPeak[++etaPeakCount] = angle;
803       }
804     }
805   } 
806
807   for (i=0; i<etaPeakCount+1; i++) {
808     fThetaCerenkov += etaPeak[i];
809   }
810   if (etaPeakCount>=0) {
811     fThetaCerenkov /= etaPeakCount+1;
812     fThetaPeakPos = etaPeakPos;
813   }
814 }
815
816
817 void AliRICHPatRec::HoughFiltering(float hcs[])
818 {
819
820 // hough filtering
821
822    float hcsFilt[750];
823    float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
824    int nx, i, nxDx;
825    int sizeHCS;
826    int nBin;
827
828    int   etaBin = 750;
829    float etaMax = 0.75;
830    float stepEta = 0.001;
831
832    nBin =  (int)(1+etaMax/stepEta); 
833    sizeHCS = etaBin*sizeof(float);
834
835    memset ((void *)hcsFilt, 0, sizeHCS); 
836
837    for (nx = 0; nx < nBin; nx++) {
838       for (i = 0; i < 5; i++)   {
839         nxDx = nx + (i-2);
840         if (nxDx> -1 && nxDx<nBin)
841              hcsFilt[nx] +=  hcs[nxDx] * k[i];
842       }      
843    }
844      
845    for (nx = 0; nx < nBin; nx++) {
846      hcs[nx] = hcsFilt[nx];
847    }
848 }
849
850 Float_t AliRICHPatRec::CherenkovRingDrawing(Float_t fixedthetacer)
851 {
852
853 //to draw Cherenkov ring by known Cherenkov angle
854
855     Int_t nmaxdegrees, nstepdegrees;
856     Float_t phpad, thetacer;
857     Float_t nfreonave, nquartzave;
858     Float_t aveEnerg;
859     Float_t energy[2];
860     Float_t e1, e2, f1, f2;
861     Float_t bandradius;
862     Float_t coordPadRing;
863     
864 //parameters to calculate freon window refractive index vs. energy
865     Float_t a = 1.177;
866     Float_t b = 0.0172;
867     
868 //parameters to calculate quartz window refractive index vs. energy
869 /*
870    Energ[0]  = 5.6;
871    Energ[1]  = 7.7;
872 */      
873     energy[0]  = 5.0;
874     energy[1]  = 8.0;
875     e1  = 10.666;
876     e2  = 18.125;
877     f1  = 46.411;
878     f2  = 228.71;
879    
880
881     nmaxdegrees = 360;
882     
883    nstepdegrees = 36;
884    
885    for (phpad=0; phpad<nmaxdegrees;phpad++) { 
886       
887        aveEnerg =  (energy[0]+energy[1])/2.;
888        
889        nfreonave  = a+b*aveEnerg;
890        nquartzave = sqrt(1+(f1/(pow(e1,2)-pow(aveEnerg,2)))+
891                          (f2/(pow(e2,2)-pow(aveEnerg,2))));
892        
893        thetacer =  fixedthetacer;
894        
895        bandradius = DistanceFromMip(nfreonave, nquartzave,
896                                    fEmissPoint,thetacer, phpad); 
897
898        coordPadRing=fPhotocatExitPhot;
899        
900        phpad = (nmaxdegrees/nstepdegrees)*phpad;
901        
902        return coordPadRing;
903    }
904
905     return coordPadRing;
906 }
907
908