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