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