Code cleaning, all wranings removed with new Makefile options
[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 /* $Id$ */
17
18
19 #include "AliRICHSDigit.h"
20 #include "AliRICHDigit.h"
21 #include "AliRICHRawCluster.h"
22 #include "AliRICHRecHit1D.h"
23 #include "AliRun.h"
24 #include "AliDetector.h"
25 #include "AliRICH.h"
26 #include "AliRICHPoints.h"
27 #include "AliSegmentation.h"
28 #include "AliRICHPatRec.h"
29 #include "AliRICH.h"
30 #include "AliRICHConst.h"
31 #include "AliRICHPoints.h"
32 #include "AliConst.h"
33 #include "AliHitMap.h"
34
35 #include <TParticle.h>
36 #include <TMath.h>
37 #include <TRandom.h>
38 #include <TCanvas.h>
39 #include <TH2.h>
40 #include <TTree.h>
41
42
43 ClassImp(AliRICHPatRec)
44 //___________________________________________
45 AliRICHPatRec::AliRICHPatRec() : TNamed()
46 {
47   // Default constructor
48   
49     //fChambers = 0;
50 }
51 //___________________________________________
52 AliRICHPatRec::AliRICHPatRec(const char *name, const char *title)
53     : TNamed(name,title)
54 {
55     //Constructor for Bari's pattern recogniton method object
56 }
57
58 void AliRICHPatRec::PatRec()
59 {
60
61 // Pattern recognition algorithm
62
63   AliRICHChamber*       iChamber;
64   AliSegmentation*  segmentation;
65         
66   Int_t ntracks, ndigits[kNCH];
67   Int_t itr, ich, i;
68   Int_t goodPhotons;
69   Int_t x,y,q;
70   Float_t rx,ry,rz;
71   Int_t nent,status;
72   Int_t padsUsedX[100];
73   Int_t padsUsedY[100];
74
75   Float_t rechit[7];
76
77   //printf("PatRec started\n");
78
79   AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
80   TTree *treeH = pRICH->TreeH();
81
82   ntracks =(Int_t) treeH->GetEntries();
83   //  ntracks = 1;
84   for (itr=0; itr<ntracks; itr++) {
85  
86     status = TrackParam(itr,ich,0,0);    
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 = &(pRICH->Chamber(ich));
92     segmentation=iChamber->GetSegmentationModel();
93
94     nent=(Int_t)gAlice->TreeD()->GetEntries();
95     gAlice->TreeD()->GetEvent(0);
96     TClonesArray *pDigitss = pRICH->DigitsAddress(ich);
97     ndigits[ich] = pDigitss->GetEntriesFast();
98     printf("Digits in chamber %d: %d\n",ich,ndigits[ich]);
99     AliRICHDigit *padI = 0;
100
101     goodPhotons = 0;
102
103     for (Int_t dig=0;dig<ndigits[ich];dig++) {
104       padI=(AliRICHDigit*) pDigitss->UncheckedAt(dig);
105       x=padI->PadX();
106       y=padI->PadY();
107       q=padI->Signal();
108       segmentation->GetPadC(x,y,rx,ry,rz);      
109
110       //printf("Pad coordinates x:%d, Real coordinates x:%f\n",x,rx);
111       //printf("Pad coordinates y:%d, Real coordinates y:%f\n",y,ry);
112
113       fXpad = rx-fXshift;
114       fYpad = ry-fYshift;
115       fQpad = q;
116     
117       fCerenkovAnglePad = PhotonCerenkovAngle();
118       if(fCerenkovAnglePad==-999) continue;
119
120       if(!PhotonInBand()) continue;
121
122       Int_t xpad;
123       Int_t ypad;
124
125       segmentation->GetPadI(fXpad,fYpad,0,xpad,ypad);
126
127       padsUsedX[goodPhotons]=xpad;
128       padsUsedY[goodPhotons]=ypad;
129       
130       goodPhotons++;
131       fEtaPhotons[goodPhotons-1] = fCerenkovAnglePad;
132     }
133     fNumEtaPhotons = goodPhotons;
134
135     BackgroundEstimation();
136
137     HoughResponse();
138     //CerenkovRingDrawing();
139
140     rechit[0] = 0;
141     rechit[1]   = 0;
142     rechit[2] = fThetaCerenkov;
143     rechit[3] = fXshift + fTrackLoc[0];
144     rechit[4] = fYshift + fTrackLoc[1];
145     rechit[5] = fEmissPoint;
146     rechit[6] = goodPhotons;
147
148     //printf("Center coordinates:%f %f\n",rechit[3],rechit[4]);
149     
150     pRICH->AddRecHit1D(ich,rechit,fEtaPhotons,padsUsedX,padsUsedY);
151     
152   }    
153
154   gAlice->TreeR()->Fill();
155   TClonesArray *fRec;
156   for (i=0;i<kNCH;i++) {
157     fRec=pRICH->RecHitsAddress1D(i);
158     int ndig=fRec->GetEntriesFast();
159     printf ("Chamber %d, rings %d\n",i,ndig);
160   }
161   pRICH->ResetRecHits1D();
162   
163 }     
164
165
166 Int_t AliRICHPatRec::TrackParam(Int_t itr, Int_t &ich, Float_t rectheta, Float_t recphi)
167 {
168   // Get Local coordinates of track impact  
169
170   AliRICHChamber*       iChamber;
171   AliSegmentation*      segmentation;
172
173   Float_t trackglob[3];
174   Float_t trackloc[3];
175   Float_t thetatr;  
176   Float_t phitr;
177   Float_t iloss;
178   Float_t part;
179   Float_t pX, pY, pZ;
180
181   //printf("Calling TrackParam\n");
182
183     gAlice->ResetHits();
184     AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
185     TTree *treeH = pRICH->TreeH();
186     treeH->GetEvent(itr);
187  
188     AliRICHhit* mHit=(AliRICHhit*)pRICH->FirstHit(-1);
189     if(mHit==0) return 1;
190     ich = mHit->Chamber()-1;
191     trackglob[0] = mHit->X();
192     trackglob[1] = mHit->Y();
193     trackglob[2] = mHit->Z();
194     pX = mHit->MomX();
195     pY = mHit->MomY();
196     pZ = mHit->MomZ();
197     fTrackMom = sqrt(TMath::Power(pX,2)+TMath::Power(pY,2)+TMath::Power(pZ,2));
198     if(recphi!=0 || rectheta!=0)
199       {
200         thetatr = rectheta;
201         phitr = recphi;
202       }
203     else
204       {
205         thetatr = mHit->Theta()*TMath::Pi()/180;
206         phitr = mHit->Phi()*TMath::Pi()/180;
207       }
208     iloss = mHit->Loss();
209     part  = mHit->Particle();
210
211     iChamber = &(pRICH->Chamber(ich));
212     iChamber->GlobaltoLocal(trackglob,trackloc);
213
214     segmentation=iChamber->GetSegmentationModel();
215
216     // retrieve geometrical params
217
218     AliRICHGeometry* fGeometry=iChamber->GetGeometryModel();   
219     
220     fRw   = fGeometry->GetFreonThickness();
221     fQw   = fGeometry->GetQuartzThickness();
222     fTgap = fGeometry->GetGapThickness(); 
223     Float_t radiatorToPads= fGeometry->GetRadiatorToPads(); 
224       //+ fGeometry->GetProximityGapThickness();
225
226     //printf("Distance to pads. From geometry:%f, From calculations:%f\n",radiatorToPads,fRw + fQw + fTgap); 
227
228     //Float_t apar = (fRw + fQw + fTgap)*tan(thetatr);
229     Float_t apar = radiatorToPads*tan(thetatr);
230     fTrackLoc[0] = apar*cos(phitr);
231     fTrackLoc[1] = apar*sin(phitr);
232     //fTrackLoc[2] = fRw + fQw + fTgap;
233     fTrackLoc[2] = radiatorToPads;
234     fTrackTheta  = thetatr;
235     fTrackPhi    = phitr;
236    
237     fXshift = trackloc[0] - fTrackLoc[0];    
238     fYshift = trackloc[2] - fTrackLoc[1];
239      
240     return 0;
241 }
242
243 Float_t AliRICHPatRec::EstimationAtLimits(Float_t lim, Float_t radius,
244                                                  Float_t phiphot)
245 {
246
247 // Estimation of emission point
248
249   Float_t nquartz = 1.585;
250   Float_t ngas    = 1.;
251   Float_t nfreon  = 1.295;
252   Float_t value;
253
254   //  printf("Calling EstimationLimits\n");
255
256   Float_t apar = (fRw -fEmissPoint + fQw + fTgap)*tan(fTrackTheta);
257   Float_t b1 = (fRw-fEmissPoint)*tan(lim);
258   Float_t b2 = fQw / sqrt(TMath::Power(nquartz,2)-TMath::Power(nfreon*sin(lim),2));
259   Float_t b3 = fTgap / sqrt(TMath::Power(ngas,2)-TMath::Power(nfreon*sin(lim),2));
260   Float_t bpar = b1 + nfreon*sin(lim)*(b2+b3);
261   value = TMath::Power(radius,2)
262     -TMath::Power((apar*cos(fTrackPhi)-bpar*cos(phiphot)),2)
263     -TMath::Power((apar*sin(fTrackPhi)-bpar*sin(phiphot)),2);
264   return value;
265 }
266
267
268 Float_t AliRICHPatRec::PhotonCerenkovAngle()
269 {
270   // Cherenkov pad angle reconstruction
271      
272   Float_t radius;
273   Float_t cherMin = 0;
274   Float_t cherMax = 0.8;
275   Float_t phiphot;
276   Float_t eps = 0.0001;
277   Int_t   niterEmiss = 0;
278   Int_t   niterEmissMax = 0;
279   Float_t x1,x2,x3=0,p1,p2,p3;
280   Float_t argY,argX;
281   Int_t niterFun;
282
283   //  printf("Calling PhotonCerenkovAngle\n");
284
285   radius = sqrt(TMath::Power(fTrackLoc[0]-fXpad,2)+TMath::Power(fTrackLoc[1]-fYpad,2));
286   fEmissPoint = fRw/2.;  //Start value of EmissionPoint
287   
288   while(niterEmiss<=niterEmissMax) {
289  
290     niterFun = 0;
291     argY = fYpad - fEmissPoint*tan(fTrackTheta)*sin(fTrackPhi);
292     argX = fXpad - fEmissPoint*tan(fTrackTheta)*cos(fTrackPhi);
293     phiphot = atan2(argY,argX); 
294     p1 = EstimationAtLimits(cherMin,radius,phiphot);
295     p2 = EstimationAtLimits(cherMax,radius,phiphot);
296     if(p1*p2>0)
297       {
298         //        printf("PhotonCerenkovAngle failed\n");
299         return -999;
300       }
301
302     //start to find the Cherenkov pad angle 
303     x1 = cherMin;
304     x2 = cherMax;
305     x3 = (x1+x2)/2.;
306     p3 = EstimationAtLimits(x3,radius,phiphot);
307     while(TMath::Abs(p3)>eps){
308       if(p1*p3<0) x2 = x3;
309       if(p1*p3>0) {
310         x1 = x3;
311         p1 = EstimationAtLimits(x1,radius,phiphot);
312       }
313       x3 = (x1+x2)/2.;
314       p3 = EstimationAtLimits(x3,radius,phiphot);
315       niterFun++;
316
317       if(niterFun>=1000) {
318         //      printf(" max iterations in PhotonCerenkovAngle\n");
319         return x3;
320       }
321     }
322     //    printf("niterFun %i \n",niterFun);
323     niterEmiss++;
324     if (niterEmiss != niterEmissMax+1) EmissionPoint();
325   } 
326   /* 
327    printf(" phiphot %f fXpad %f fYpad %f fEmiss %f \n",
328                          phiphot,fXpad,fYpad,fEmissPoint);
329   */
330
331   return x3;
332
333 }
334
335
336 void AliRICHPatRec::EmissionPoint()
337
338
339 // Find emission point
340
341   Float_t absorbtionLength=7.83*fRw; //absorption length in the freon (cm)
342   // 7.83 = -1/ln(T0) where 
343   // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
344   Float_t photonLength, photonLengthMin, photonLengthMax;
345
346   photonLength=exp(-fRw/(absorbtionLength*cos(fCerenkovAnglePad)));
347   photonLengthMin=fRw*photonLength/(1.-photonLength);
348   photonLengthMax=absorbtionLength*cos(fCerenkovAnglePad);
349   fEmissPoint = fRw + photonLengthMin - photonLengthMax;
350
351 }
352
353 void AliRICHPatRec::PhotonSelection(Int_t track, Int_t &nphot, Float_t &thetamean)
354 {
355
356 // not implemented yet
357   track++;
358   nphot=10;
359   thetamean=4;
360   printf("Calling PhotonSelection\n");
361 }
362
363 void AliRICHPatRec::BackgroundEstimation()
364 {
365
366 // estimate background noise
367
368   Float_t stepEta   = 0.001;  
369   Float_t etaMinBkg = 0.72;
370   Float_t etaMaxBkg = 0.75;
371   Float_t etaMin    = 0.;
372   Float_t etaMax    = 0.75;
373   Float_t ngas      = 1.;
374   Float_t nfreon    = 1.295;
375
376   Float_t etaStepMin,etaStepMax,etaStepAvg;
377   Int_t i,ip,nstep;
378   Int_t numPhotBkg, numPhotonStep;
379   Float_t funBkg,areaBkg,normBkg;
380   Float_t densityBkg,storeBkg,numStore;
381   Float_t thetaSig;
382   
383   numPhotBkg = 0;
384   areaBkg = 0.;
385
386   nstep = (int)((etaMaxBkg-etaMinBkg)/stepEta);
387
388   for (i=0;i<fNumEtaPhotons;i++) {
389
390     if(fEtaPhotons[i]>etaMinBkg && fEtaPhotons[i]<etaMaxBkg) {
391       numPhotBkg++;
392     }    
393   }
394   if (numPhotBkg == 0) {
395      for (i=0;i<fNumEtaPhotons;i++) {
396         fWeightPhotons[i] = 1.;
397      }
398     return;
399   }
400
401   //  printf(" numPhotBkg %i ",numPhotBkg);
402
403   for (i=0;i<nstep;i++) {
404     etaStepMin = etaMinBkg + (Float_t)(i)*stepEta;
405     etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;    
406     etaStepAvg = 0.5*(etaStepMax + etaStepMin);
407     /*
408     funBkg = tan(etaStepAvg)*TMath::Power((1.+TMath::Power(tan(etaStepAvg),2)),
409                                   5.52)-7.803 + 22.02*tan(etaStepAvg);
410     */
411     
412     //printf("etaStepAvg: %f, etaStepMax: %f, etaStepMin: %f", etaStepAvg,etaStepMax,etaStepMin);
413
414     thetaSig = TMath::ASin(nfreon/ngas*TMath::Sin(etaStepAvg));
415     funBkg = tan(thetaSig)*(1.+TMath::Power(tan(thetaSig),2))*nfreon
416        /ngas*cos(etaStepAvg)/cos(thetaSig);
417     areaBkg += stepEta*funBkg;
418   }
419
420   densityBkg = 0.95*(Float_t)(numPhotBkg)/areaBkg;
421   //  printf(" densityBkg %f \n",densityBkg);
422   
423   nstep = (int)((etaMax-etaMin)/stepEta); 
424   storeBkg = 0.;
425   numStore = 0;
426   for (i=0;i<nstep;i++) {
427     etaStepMin = etaMinBkg + (Float_t)(i)*stepEta;
428     etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;    
429     etaStepAvg = 0.5*(etaStepMax + etaStepMin);
430     /*
431     funBkg = tan(etaStepAvg)*TMath::Power((1.+TMath::Power(tan(etaStepAvg),2)),
432                                   5.52)-7.803 + 22.02*tan(etaStepAvg);
433     */
434
435     thetaSig = asin(nfreon/ngas*sin(etaStepAvg));
436     funBkg = tan(thetaSig)*(1.+TMath::Power(tan(thetaSig),2))*nfreon
437        /ngas*cos(etaStepAvg)/cos(thetaSig);
438
439     areaBkg = stepEta*funBkg;
440     normBkg = densityBkg*areaBkg; 
441     numPhotonStep = 0;
442     for (ip=0;ip<fNumEtaPhotons;ip++) {
443       if(fEtaPhotons[ip]>etaStepMin && fEtaPhotons[ip]<etaStepMax) {
444         numPhotonStep++;
445       }
446     }
447     if (numPhotonStep == 0) {
448       storeBkg += normBkg;
449       numStore++;
450       if (numStore>50) {
451         numStore = 0;
452         storeBkg = 0.;
453       }
454     }
455     if (numPhotonStep == 0) continue;
456     for (ip=0;ip<fNumEtaPhotons;ip++) {
457       if(fEtaPhotons[ip]>etaStepMin && fEtaPhotons[ip]<etaStepMax) {
458         normBkg +=storeBkg;
459         storeBkg = 0;
460         numStore = 0;
461         fWeightPhotons[ip] = 1. - normBkg/(Float_t)(numPhotonStep);
462         /*
463         printf(" normBkg %f numPhotonStep %i fW %f \n",
464                normBkg, numPhotonStep, fWeightPhotons[ip]);
465         */
466         if(fWeightPhotons[ip]<0) fWeightPhotons[ip] = 0.;
467       }
468     }
469   }
470 }
471
472
473 void AliRICHPatRec::FlagPhotons(Int_t track, Float_t theta)
474 {
475
476 // not implemented yet
477   track++;theta++;
478   printf("Calling FlagPhotons\n");
479 }
480
481
482 //////////////////////////////////////////
483
484
485
486
487
488 Int_t AliRICHPatRec::PhotonInBand()
489
490   //0=label for parameters giving internal band ellipse
491   //1=label for parameters giving external band ellipse  
492
493   Float_t imp[2], mass[2], energy[2], beta[2]; 
494   Float_t emissPointLength[2];
495   Float_t e1, e2, f1, f2;
496   Float_t nfreon[2], nquartz[2]; 
497   Int_t times;
498   Float_t pointsOnCathode[3];
499
500   Float_t phpad, thetacer[2]; 
501   Float_t bandradius[2], padradius;
502
503   imp[0] = 5.0; //threshold momentum for the proton Cherenkov emission
504   imp[1] = 1.2;
505  
506   mass[0] = 0.938; //proton mass 
507   mass[1] = 0.139; //pion mass
508
509   emissPointLength[0] = fRw-0.0001; //at the beginning of the radiator
510   emissPointLength[1] = 0.;//at the end of radiator
511   
512   //parameters to calculate freon window refractive index vs. energy
513   Float_t a = 1.177;
514   Float_t b = 0.0172;
515  
516   //parameters to calculate quartz window refractive index vs. energy
517   /*
518   Energ[0]  = 5.6;
519   Energ[1]  = 7.7;
520   */
521   energy[0]  = 5.0;
522   energy[1]  = 8.0;
523   e1 = 10.666;
524   e2  = 18.125;
525   f1  = 46.411;
526   f2  = 228.71;
527
528   phpad = PhiPad(fTrackTheta,fTrackPhi);  
529
530   for (times=0; times<=1; times++) {
531   
532     nfreon[times]   = a+b*energy[times];
533     //nfreon[times]   = 1;
534
535     nquartz[times] = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(energy[times],2)))+
536                           (f2/(TMath::Power(e2,2)-TMath::Power(energy[times],2))));
537
538     beta[times]  = imp[times]/sqrt(TMath::Power(imp[times],2)+TMath::Power(mass[times],2));
539    
540     thetacer[times] =  CherenkovAngle( nfreon[times], beta[times]);
541
542     bandradius[times] = DistanceFromMip( nfreon[times], nquartz[times],
543                                 emissPointLength[times], 
544                                       thetacer[times], phpad, pointsOnCathode,fTrackTheta,fTrackPhi);
545     //printf(" ppp %f %f %f \n",pointsOnCathode);  
546 }
547
548   bandradius[0] -= 1.6;
549   bandradius[1] += 1.6;
550   padradius = sqrt(TMath::Power(fXpad,2)+TMath::Power(fYpad,2));
551   //  printf(" rmin %f r %f rmax %f \n",bandradius[0],padradius,bandradius[1]);
552
553   if(padradius>=bandradius[0] && padradius<=bandradius[1]) return 1;
554   return 0; 
555 }
556
557 Float_t AliRICHPatRec::DistanceFromMip(Float_t nfreon, Float_t nquartz, 
558                        Float_t emissPointLength, Float_t thetacer, 
559                        Float_t phpad, Float_t pointsOnCathode[3], Float_t rectheta, Float_t recphi)
560
561
562 // Find the distance to MIP impact
563
564   Float_t distanceValue;
565
566   TVector3 radExitPhot(1,1,1);//photon impact at the radiator exit with respect
567   //to local reference sistem with the origin in the MIP entrance 
568    
569   TVector3 vectEmissPointLength(1,1,1);
570   Float_t magEmissPointLenght;
571
572   TVector3 radExitPhot2(1,1,1);//photon impact at the radiator exit with respect
573   Float_t magRadExitPhot2;
574   //to a reference sistem with origin in the photon emission point and  
575   //axes parallel to the MIP reference sistem
576
577   TVector3 quarExitPhot(1,1,1);//photon impact at the quartz exit with respect
578   Float_t magQuarExitPhot;
579   // 
580   TVector3 gapExitPhot(1,1,1) ;
581   Float_t magGapExitPhot;
582   //
583   TVector3 PhotocatExitPhot(1,1,1);
584   Double_t theta2;
585   Double_t thetarad , phirad ;
586   Double_t thetaquar, phiquar;
587   Double_t thetagap , phigap ;
588
589   Float_t ngas    = 1.;
590
591   magEmissPointLenght =  emissPointLength/cos(rectheta);
592
593   vectEmissPointLength.SetMag(magEmissPointLenght);
594   vectEmissPointLength.SetTheta(rectheta);
595   vectEmissPointLength.SetPhi(recphi);
596
597
598   radExitPhot2.SetTheta(thetacer);  
599   radExitPhot2.SetPhi(phpad); 
600
601
602   TRotation r1;
603   TRotation r2;
604   TRotation r;
605
606   r1. RotateY(rectheta);
607   r2. RotateZ(recphi);
608   
609
610
611   r = r2 * r1;//rotation about the y axis by MIP theta incidence angle
612   //following by a rotation about the z axis by MIP phi incidence angle; 
613
614
615   radExitPhot2    = r * radExitPhot2;
616   theta2          = radExitPhot2.Theta();
617   magRadExitPhot2 = (fRw -  vectEmissPointLength(2))/cos(theta2);
618   radExitPhot2.SetMag(magRadExitPhot2);
619
620
621   radExitPhot = vectEmissPointLength + radExitPhot2;
622   thetarad    = radExitPhot.Theta();
623   phirad  =  radExitPhot.Phi(); //check on the original file //
624
625   thetaquar   = SnellAngle( nfreon, nquartz, theta2); 
626   phiquar     = radExitPhot2.Phi(); 
627   if(thetaquar == 999.) return thetaquar;
628   magQuarExitPhot    = fQw/cos(thetaquar);
629   quarExitPhot.SetMag( magQuarExitPhot);
630   quarExitPhot.SetTheta(thetaquar);
631   quarExitPhot.SetPhi(phiquar);
632
633   thetagap = SnellAngle( nquartz, ngas, thetaquar); 
634   phigap   = phiquar; 
635   if(thetagap == 999.) return thetagap;
636   magGapExitPhot    = fTgap/cos(thetagap);
637   gapExitPhot.SetMag( magGapExitPhot);
638   gapExitPhot.SetTheta(thetagap);
639   gapExitPhot.SetPhi(phigap);
640
641   PhotocatExitPhot =  radExitPhot + quarExitPhot + gapExitPhot; 
642
643   distanceValue = sqrt(TMath::Power(PhotocatExitPhot(0),2)
644                            +TMath::Power(PhotocatExitPhot(1),2)); 
645   pointsOnCathode[0] = (Float_t) PhotocatExitPhot(0) + fXshift - fTrackLoc[0];
646   pointsOnCathode[1] = (Float_t) PhotocatExitPhot(1) + fYshift - fTrackLoc[1];
647   pointsOnCathode[2] = (Float_t) PhotocatExitPhot(2);
648
649   //printf(" point in Distance.2. %f %f %f \n",pointsOnCathode[0],pointsOnCathode[1],pointsOnCathode[2]);
650
651   return distanceValue;
652
653  }
654
655 Float_t AliRICHPatRec::PhiPad(Float_t rectheta, Float_t recphi)
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 = - rectheta;
675   phirot   = - recphi;
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 /*void AliRICHPatRec::CerenkovRingDrawing()
851 {
852
853 //to draw Cherenkov ring by known Cherenkov angle
854
855     Int_t nmaxdegrees;
856     Int_t Nphpad;
857     Float_t phpad;
858     Float_t nfreonave, nquartzave;
859     Float_t aveEnerg;
860     Float_t energy[2];
861     Float_t e1, e2, f1, f2;
862     Float_t bandradius;
863
864 //parameters to calculate freon window refractive index vs. energy
865
866     Float_t a = 1.177;
867     Float_t b = 0.0172;
868     
869 //parameters to calculate quartz window refractive index vs. energy
870
871     energy[0]  = 5.0;
872     energy[1]  = 8.0;
873     e1  = 10.666;
874     e2  = 18.125;
875     f1  = 46.411;
876     f2  = 228.71;
877    
878
879    nmaxdegrees = 36;
880    
881    for (Nphpad=0; Nphpad<nmaxdegrees;Nphpad++) { 
882
883        phpad = (360./(Float_t)nmaxdegrees)*(Float_t)Nphpad;
884       
885        aveEnerg =  (energy[0]+energy[1])/2.;
886        
887        nfreonave  = a+b*aveEnerg;
888        nquartzave = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(aveEnerg,2)))+
889                          (f2/(TMath::Power(e2,2)-TMath::Power(aveEnerg,2))));
890        
891        bandradius = DistanceFromMip(nfreonave, nquartzave,
892                                    fEmissPoint,fThetaCerenkov, phpad);
893
894        fCoordEllipse[0][Nphpad] = fOnCathode[0];
895        fCoordEllipse[1][Nphpad] = fOnCathode[1];
896        printf(" values %f %f \n",fOnCathode[0],fOnCathode[1]);
897        
898    }
899
900 }*/
901
902