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