b3021fda02972c627145c92a2e34d158dbd3486d
[u/mrichter/AliRoot.git] / RICH / AliRICHRecon.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 // AliRICHRecon                                                         //
19 //                                                                      //
20 // RICH class to perfom pattern recognition based on Hough transfrom    //
21 //                                                                      //
22 //////////////////////////////////////////////////////////////////////////
23
24 #include <Riostream.h>
25 #include <TMath.h>
26 #include <TRotation.h>
27 #include <TVector3.h>
28
29 #include "AliRICH.h"
30 #include "AliRICHParam.h"
31 #include "AliRICHRecon.h"
32 #include "AliRICHHelix.h"
33 #include <AliLog.h>
34
35 #define NPointsOfRing 201
36
37 //__________________________________________________________________________________________________
38 AliRICHRecon::AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId)
39              :TTask("RichRec","RichPat")
40 {
41 // main ctor
42   SetFreonScaleFactor(1);
43   fIsWEIGHT = kFALSE;
44   fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75; 
45   fDTheta       = 0.001;   fWindowWidth  = 0.045;
46   fRadiatorWidth = AliRICHParam::Zfreon();
47   fQuartzWidth   = AliRICHParam::Zwin();
48   fGapWidth      = AliRICHParam::Freon2Pc() - fRadiatorWidth - fQuartzWidth;
49   fXmin         = -AliRICHParam::PcSizeX()/2.;
50   fXmax         =  AliRICHParam::PcSizeX()/2.;
51   fYmin         = -AliRICHParam::PcSizeY()/2.;
52   fYmax         =  AliRICHParam::PcSizeY()/2.; 
53   SetTrackTheta(pHelix->Ploc().Theta());
54   SetTrackPhi(pHelix->Ploc().Phi());
55   SetMipIndex(iMipId);
56   SetShiftX(pHelix->PosRad().X());
57   SetShiftY(pHelix->PosRad().Y());
58   fpClusters = pClusters;
59 }
60 //__________________________________________________________________________________________________
61 Double_t AliRICHRecon::ThetaCerenkov()
62 {
63 // Pattern recognition method based on Hough transform
64 // Return theta Cerenkov for a given track and list of clusters which are set in ctor  
65
66   if(fpClusters->GetEntries()==0) return -1;//no clusters at all for a given track
67   Bool_t kPatRec = kFALSE;  
68     
69   AliDebug(1,Form("---Track Parameters--- Theta: %f , Phi: %f ",GetTrackTheta()*TMath::RadToDeg(),GetTrackPhi()*TMath::RadToDeg()));
70
71   Int_t candidatePhotons = 0;
72
73   SetThetaCerenkov(999.);
74   SetHoughPhotons(0);
75   SetHoughPhotonsNorm(0);
76
77   for (Int_t j=0; j < fpClusters->GetEntries(); j++){//clusters loop
78     SetPhotonIndex(j);
79     SetPhotonFlag(0);
80     SetPhotonEta(-999.);
81     SetPhotonWeight(0.);
82     if (j == GetMipIndex()) continue; // do not consider MIP cluster as a candidate photon
83     Float_t xtoentr = ((AliRICHCluster*)fpClusters->UncheckedAt(j))->X() - GetShiftX();
84     Float_t ytoentr = ((AliRICHCluster*)fpClusters->UncheckedAt(j))->Y() - GetShiftY();
85     SetEntranceX(xtoentr);
86     SetEntranceY(ytoentr);
87     FindPhiPoint();
88 //      Int_t photonStatus = PhotonInBand();
89 //      if(photonStatus == 0) continue;
90     SetPhotonFlag(1);
91     FindThetaPhotonCerenkov();
92     Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
93     AliDebug(1,Form("THETA CERENKOV ---> %f",thetaPhotonCerenkov));
94     SetPhotonEta(thetaPhotonCerenkov);
95     candidatePhotons++;
96   }//clusters loop
97
98   if(candidatePhotons >= 1) kPatRec = kTRUE;
99
100   if(!kPatRec) return -1;
101
102   SetPhotonsNumber(fpClusters->GetEntries());
103
104   HoughResponse();
105   
106   fNrings++;
107
108   FlagPhotons();
109   Int_t nPhotonHough = GetHoughPhotons();
110  
111   if(nPhotonHough < 1) 
112     {
113       SetThetaCerenkov(999.);
114       SetHoughPhotonsNorm(0.);
115       return -1;
116     }
117
118   FindThetaCerenkov();
119
120   AliDebug(1,Form("Number of clusters accepted --->  %i",nPhotonHough));
121   
122 //  Float_t thetaCerenkov = GetThetaCerenkov();  
123 //  SetThetaOfRing(thetaCerenkov);
124 //  FindAreaAndPortionOfRing();
125
126 //  Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
127 //  SetHoughPhotonsNorm(nPhotonHoughNorm);
128
129   // Calculate the area where the photon are accepted...
130 /*
131   Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth; 
132   SetThetaOfRing(thetaInternal);
133   FindAreaAndPortionOfRing();
134   Float_t internalArea = GetAreaOfRing();
135
136   Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth; 
137   SetThetaOfRing(thetaExternal);
138   FindAreaAndPortionOfRing();
139   Float_t externalArea = GetAreaOfRing();
140
141   Float_t houghArea = externalArea - internalArea;
142
143   SetHoughArea(houghArea);
144 */
145   return GetThetaCerenkov();
146
147 }//ThetaCerenkov()
148 //__________________________________________________________________________________________________
149 void AliRICHRecon::FindEmissionPoint()
150 {
151   //estimate the emission point in radiator
152
153 // Find emission point
154
155   Float_t absorbtionLenght=7.83*fRadiatorWidth; //absorption length in the freon (cm)
156   // 7.83 = -1/ln(T0) where 
157   // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
158   Float_t photonLenght, photonLenghtMin, photonLenghtMax;
159
160   photonLenght=exp(-fRadiatorWidth/(absorbtionLenght*cos(fCerenkovAnglePad)));
161   photonLenghtMin=fRadiatorWidth*photonLenght/(1.-photonLenght);
162   photonLenghtMax=absorbtionLenght*cos(fCerenkovAnglePad);
163   Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax;
164
165   SetEmissionPoint(emissionPoint);
166   SetEmissionPoint(fRadiatorWidth/2); // tune the emission point
167 }
168 //__________________________________________________________________________________________________
169 Int_t AliRICHRecon::PhotonInBand()
170 {
171   //search band fro photon candidates
172
173   //  Float_t massOfParticle;
174   Float_t nfreon;
175
176   Float_t thetacer;
177
178   Float_t xtoentr = GetEntranceX();
179   Float_t ytoentr = GetEntranceY();
180
181   Float_t innerRadius;
182   Float_t outerRadius;
183
184   Float_t phpad = GetPhiPoint();
185
186
187   // inner radius //
188   SetPhotonEnergy(5.6);
189   SetEmissionPoint(fRadiatorWidth -0.0001);
190   SetFreonRefractiveIndex();
191
192   nfreon = GetFreonRefractiveIndex();
193   thetacer = 0.;
194
195   AliDebug(1,Form("thetacer in photoninband min %f",thetacer));
196
197   FindThetaAtQuartz(thetacer);
198
199   if(thetacer == 999. || GetThetaAtQuartz() == 999.)
200     {
201       innerRadius = -999.;
202       SetXInnerRing(-999.);
203       SetYInnerRing(-999.);
204       SetRadiusInnerRing(-999.);
205     }
206   else
207     {
208       SetThetaPhotonInDRS(GetThetaAtQuartz());
209       SetPhiPhotonInDRS(phpad);
210
211       innerRadius = FromEmissionToCathode();
212        if(innerRadius == 999.) innerRadius = -999.;
213       
214       SetXInnerRing(GetXPointOnCathode());
215       SetYInnerRing(GetYPointOnCathode());
216       SetRadiusInnerRing(innerRadius);
217     }
218   
219   // outer radius //
220   SetPhotonEnergy(7.7);
221   SetEmissionPoint(0.);
222 //  SetMassHypotesis(0.139567);
223   SetFreonRefractiveIndex();
224
225   nfreon = GetFreonRefractiveIndex();
226
227   thetacer = Cerenkovangle(nfreon,1);
228
229   //  thetacer = 0.75;
230
231   AliDebug(1,Form("thetacer in photoninband max %f",thetacer));
232
233   FindThetaAtQuartz(thetacer);
234
235   if(thetacer == 999. || GetThetaAtQuartz() == 999.)
236     {
237       outerRadius = 999.;
238       SetXOuterRing(999.);
239       SetYOuterRing(999.);
240       SetRadiusOuterRing(999.);
241     }
242   else
243     {
244       SetThetaPhotonInDRS(GetThetaAtQuartz());
245       SetPhiPhotonInDRS(phpad);
246
247       outerRadius = FromEmissionToCathode();
248 //      cout << " outerRadius " << outerRadius << endl;
249       SetXOuterRing(GetXPointOnCathode());
250       SetYOuterRing(GetYPointOnCathode());
251       SetRadiusOuterRing(outerRadius);
252     }
253
254   Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2));
255   
256   AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius));
257
258   if(padradius>=innerRadius && padradius<=outerRadius) return 1;
259   return 0;
260 }
261 //__________________________________________________________________________________________________
262 void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
263 {
264   //find the theta at the quartz plate
265
266   if(thetaCerenkov == 999.) 
267     {
268       SetThetaAtQuartz(999.);
269       return;
270     }
271
272   Float_t thetaAtQuartz = 999.;
273
274   Float_t trackTheta = GetTrackTheta();
275
276   if(trackTheta == 0) {
277
278     thetaAtQuartz = thetaCerenkov;
279     SetThetaAtQuartz(thetaAtQuartz);
280     return;
281   }
282
283   Float_t trackPhi   = GetTrackPhi();
284   Float_t phiPoint = GetPhiPoint();
285
286   Double_t den = TMath::Sin((Double_t)trackTheta)
287     *TMath::Cos((Double_t)trackPhi)
288     *TMath::Cos((Double_t)phiPoint) +
289     TMath::Sin((Double_t)trackTheta)
290     *TMath::Sin((Double_t)trackPhi)
291     *TMath::Sin((Double_t)phiPoint); 
292   Double_t b = TMath::Cos((Double_t)trackTheta)/den;
293   Double_t c = -TMath::Cos((Double_t)thetaCerenkov)/den;
294
295   Double_t underSqrt = 1 + b*b - c*c;
296
297   if(underSqrt < 0) {
298     SetThetaAtQuartz(999.);
299     return;
300   }
301
302   Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c);
303   Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c);
304
305   Double_t thetaSol1 = 2*TMath::ATan(sol1);
306   Double_t thetaSol2 = 2*TMath::ATan(sol2);
307
308   if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1;
309   if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2;
310
311 //  AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz));
312
313   SetThetaAtQuartz(thetaAtQuartz);
314 }
315 //__________________________________________________________________________________________________
316 void AliRICHRecon::FindThetaPhotonCerenkov()
317 {
318   //find theta cerenkov of ring
319
320   Float_t thetaCerMin = 0.;
321   Float_t thetaCerMax = 0.75;
322   Float_t thetaCerMean;
323
324   Float_t radiusMin, radiusMax, radiusMean;
325   Int_t nIteration = 0;
326
327   const Float_t kTollerance = 0.05;
328
329
330   Float_t phiPoint = GetPhiPoint();
331
332   SetPhotonEnergy(AliRICHParam::MeanCkovEnergy());
333   SetEmissionPoint(fRadiatorWidth/2);
334
335   Float_t xPoint = GetEntranceX();
336   Float_t yPoint = GetEntranceY();
337   Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint);
338
339 //  AliDebug(1,Form(" DistPoint %f ",distPoint));
340
341   // Star minimization...
342
343   // First value...
344
345   FindThetaAtQuartz(thetaCerMin);
346   
347   if(GetThetaAtQuartz() == 999.)
348     {
349       radiusMin = -999.;
350     }
351   else
352     {
353       SetThetaPhotonInDRS(GetThetaAtQuartz());
354       SetPhiPhotonInDRS(phiPoint);
355       
356       radiusMin = FromEmissionToCathode();
357     }
358
359   // Second value...
360
361   FindThetaAtQuartz(thetaCerMax);
362   if(GetThetaAtQuartz() == 999.)
363     {
364       radiusMax = 999.;
365     }
366   else
367     {
368       SetThetaPhotonInDRS(GetThetaAtQuartz());
369       SetPhiPhotonInDRS(phiPoint);
370       
371       radiusMax = FromEmissionToCathode();
372     }
373   // Mean value...
374
375   thetaCerMean = (thetaCerMax + thetaCerMin)/2;
376
377   FindThetaAtQuartz(thetaCerMean);
378   if(GetThetaAtQuartz() == 999.)
379     {
380       radiusMean = 999.;
381     }
382   else
383     {
384       SetThetaPhotonInDRS(GetThetaAtQuartz());
385       SetPhiPhotonInDRS(phiPoint);
386       
387       radiusMean = FromEmissionToCathode();
388     }
389
390 //  AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax));
391
392   while (TMath::Abs(radiusMean-distPoint) > kTollerance)
393     {
394
395       if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean;
396       if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) {
397
398         thetaCerMin = thetaCerMean;
399
400         FindThetaAtQuartz(thetaCerMin);
401         SetThetaPhotonInDRS(GetThetaAtQuartz());
402         SetPhiPhotonInDRS(phiPoint);
403
404         radiusMin =FromEmissionToCathode();
405       }
406
407       thetaCerMean = (thetaCerMax + thetaCerMin)/2;
408
409       FindThetaAtQuartz(thetaCerMean);
410       SetThetaPhotonInDRS(GetThetaAtQuartz());
411       SetPhiPhotonInDRS(phiPoint);
412
413       radiusMean = FromEmissionToCathode();
414
415       nIteration++;
416       if(nIteration>=50) {
417 //      AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration));
418         SetThetaPhotonCerenkov(999.);
419         return;
420       }
421     }
422
423 //  AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean));
424   SetThetaPhotonCerenkov(thetaCerMean);
425
426 }
427 //__________________________________________________________________________________________________
428 void AliRICHRecon::FindAreaAndPortionOfRing()
429 {
430   //find fraction of the ring accepted by the RICH
431
432   Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing];
433
434   //  Float_t xtoentr = GetEntranceX();
435   //  Float_t ytoentr = GetEntranceY();
436   Float_t shiftX = GetShiftX();
437   Float_t shiftY = GetShiftY();
438
439   Float_t xemiss = GetXCoordOfEmission(); 
440   Float_t yemiss = GetYCoordOfEmission(); 
441
442   Float_t x0 = xemiss + shiftX;
443   Float_t y0 = yemiss + shiftY;
444
445
446   SetPhotonEnergy(AliRICHParam::MeanCkovEnergy());
447   SetFreonRefractiveIndex();
448
449   SetEmissionPoint(fRadiatorWidth/2.);
450
451   Float_t theta = GetThetaOfRing();
452   
453   Int_t nPoints = 0;
454   Int_t nPsiAccepted = 0;
455   Int_t nPsiTotal = 0;
456
457   for(Int_t i=0;i<NPointsOfRing-1;i++)
458     {
459
460       Float_t psi = 2*TMath::Pi()*i/NPointsOfRing;
461       
462       SetThetaPhotonInTRS(theta);
463       SetPhiPhotonInTRS(psi);
464       FindPhotonAnglesInDRS();
465       
466       Float_t radius = FromEmissionToCathode();
467       if (radius == 999.) continue;
468       
469       nPsiTotal++;
470
471       Float_t xPointRing = GetXPointOnCathode() + shiftX;
472       Float_t yPointRing = GetYPointOnCathode() + shiftY;
473       
474       SetDetectorWhereX(xPointRing);
475       SetDetectorWhereY(yPointRing);
476       
477       Int_t zone = CheckDetectorAcceptance();
478
479
480       if (zone != 0) 
481         {
482           FindIntersectionWithDetector();
483           xPoint[nPoints] = GetIntersectionX();
484           yPoint[nPoints] = GetIntersectionY();
485         }
486       else
487         {
488           xPoint[nPoints] = xPointRing;
489           yPoint[nPoints] = yPointRing;
490           nPsiAccepted++;
491         }
492
493       nPoints++;
494
495     }
496
497   xPoint[nPoints] = xPoint[0];
498   yPoint[nPoints] = yPoint[0];
499   
500   // find area...
501
502   Float_t area = 0;
503
504   for (Int_t i = 0; i < nPoints; i++)
505     {
506       area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0));
507     }
508   
509   area *= 0.5;
510   
511   Float_t portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal));
512
513
514   SetAreaOfRing(area);
515   SetPortionOfRing(portionOfRing);
516 }
517 //__________________________________________________________________________________________________
518 void AliRICHRecon::FindIntersectionWithDetector()
519 {
520   // find ring intersection with CsI edges
521
522   Float_t xIntersect, yIntersect;
523   Float_t x1, x2, y1, y2;
524
525   Float_t shiftX = GetShiftX();
526   Float_t shiftY = GetShiftY();
527
528   Float_t xPoint = GetXPointOnCathode() + shiftX;
529   Float_t yPoint = GetYPointOnCathode() + shiftY;
530
531   Float_t xemiss = GetXCoordOfEmission(); 
532   Float_t yemiss = GetYCoordOfEmission(); 
533
534   Float_t phi = GetPhiPhotonInDRS();
535   Float_t m = tan(phi);
536
537   Float_t x0 = xemiss + shiftX;
538   Float_t y0 = yemiss + shiftY;
539
540   if(xPoint > x0)
541     {
542       x1 = x0;
543       x2 = xPoint;
544     }
545   else
546     {
547       x2 = x0;
548       x1 = xPoint;
549     }
550   if(yPoint > y0)
551     {
552       y1 = y0;
553       y2 = yPoint;
554     }
555   else
556     {
557       y2 = y0;
558       y1 = yPoint;
559     }
560   //
561   xIntersect = fXmax;
562   yIntersect = m*(xIntersect - x0) + y0;
563   if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
564     {
565       SetIntersectionX(xIntersect);
566       SetIntersectionY(yIntersect);
567       return;
568     }
569   //
570   xIntersect = fXmin;
571   yIntersect = m*(xIntersect - x0) + y0;
572   if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
573     {
574       SetIntersectionX(xIntersect);
575       SetIntersectionY(yIntersect);
576       return;
577     }
578   //
579   yIntersect = fYmax;
580   xIntersect = (yIntersect - y0)/m + x0;
581   if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
582     {
583       SetIntersectionX(xIntersect);
584       SetIntersectionY(yIntersect);
585       return;
586     }
587   //
588   yIntersect = fYmin;
589   xIntersect = (yIntersect - y0)/m + x0;
590   if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
591     {
592       SetIntersectionX(xIntersect);
593       SetIntersectionY(yIntersect);
594       return;
595     }
596   
597   cout << " sono fuori!!!!!!" << endl;
598   
599 }
600 //__________________________________________________________________________________________________
601 Int_t AliRICHRecon::CheckDetectorAcceptance() const
602 {
603   // check for the acceptance
604
605   // crosses X -2.6 2.6 cm
606   // crosses Y -1 1 cm
607
608   Float_t xcoord = GetDetectorWhereX();
609   Float_t ycoord = GetDetectorWhereY();
610
611   if(xcoord > fXmax)
612     {
613       if(ycoord > fYmax) return 2;
614       if(ycoord > fYmin && ycoord < fYmax) return 3;
615       if(ycoord < fYmin) return 4;
616     }
617   if(xcoord < fXmin)
618     {
619       if(ycoord > fYmax) return 8;
620       if(ycoord > fYmin && ycoord < fYmax) return 7;
621       if(ycoord < fYmin) return 6;
622     }
623   if(xcoord > fXmin && xcoord < fXmax)
624     {
625       if(ycoord > fYmax) return 1;
626       if(ycoord > fYmin && ycoord < fYmax) return 0;
627       if(ycoord < fYmin) return 5;
628     }
629   return 999;
630 }
631 //__________________________________________________________________________________________________
632 void AliRICHRecon::FindPhotonAnglesInDRS()
633 {
634   // Setup the rotation matrix of the track...
635
636   TRotation mtheta;
637   TRotation mphi;
638   TRotation minv;
639   TRotation mrot;
640   
641   Float_t trackTheta = GetTrackTheta();
642   Float_t trackPhi = GetTrackPhi();
643
644   mtheta.RotateY(trackTheta);
645   mphi.RotateZ(trackPhi);
646   
647   mrot = mphi * mtheta;
648   //  minv = mrot.Inverse();
649
650   TVector3 photonInRadiator(1,1,1);
651
652   Float_t thetaCerenkov = GetThetaPhotonInTRS();
653   Float_t phiCerenkov   = GetPhiPhotonInTRS();
654
655   photonInRadiator.SetTheta(thetaCerenkov);
656   photonInRadiator.SetPhi(phiCerenkov);
657   photonInRadiator = mrot * photonInRadiator;
658   Float_t theta = photonInRadiator.Theta();
659   Float_t phi = photonInRadiator.Phi();
660   SetThetaPhotonInDRS(theta);
661   SetPhiPhotonInDRS(phi);
662
663 }
664 //__________________________________________________________________________________________________
665 Float_t AliRICHRecon::FromEmissionToCathode()
666 {
667   // trace from emission point to cathode
668
669   Float_t nfreon, nquartz, ngas; 
670
671   SetFreonRefractiveIndex();
672   SetQuartzRefractiveIndex();
673   SetGasRefractiveIndex();
674
675   nfreon  = GetFreonRefractiveIndex();
676   nquartz = GetQuartzRefractiveIndex();
677   ngas    = GetGasRefractiveIndex();
678
679   Float_t trackTheta = GetTrackTheta();
680   Float_t trackPhi = GetTrackPhi();
681   Float_t lengthOfEmissionPoint = GetEmissionPoint();
682
683   Float_t theta = GetThetaPhotonInDRS();
684   Float_t phi   = GetPhiPhotonInDRS();
685
686 //   cout << " Theta " << Theta << " Phi " << Phi << endl;
687
688   Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi);
689   Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi);
690
691   SetXCoordOfEmission(xemiss);
692   SetYCoordOfEmission(yemiss);
693   
694   Float_t thetaquar = SnellAngle(nfreon, nquartz, theta);
695
696   if(thetaquar == 999.) 
697     {
698       SetXPointOnCathode(999.);
699       SetYPointOnCathode(999.);
700       return thetaquar;
701     }
702
703   Float_t thetagap  = SnellAngle( nquartz, ngas, thetaquar);
704
705   if(thetagap == 999.) 
706     {
707       SetXPointOnCathode(999.);
708       SetYPointOnCathode(999.);
709       return thetagap;
710     }
711
712   Float_t xw = (fRadiatorWidth - lengthOfEmissionPoint)*cos(phi)*tan(theta);
713   Float_t xq = fQuartzWidth*cos(phi)*tan(thetaquar);
714   Float_t xg = fGapWidth*cos(phi)*tan(thetagap);
715   Float_t yw = (fRadiatorWidth - lengthOfEmissionPoint)*sin(phi)*tan(theta);
716   Float_t yq = fQuartzWidth*sin(phi)*tan(thetaquar);
717   Float_t yg = fGapWidth*sin(phi)*tan(thetagap);
718
719
720   Float_t xtot = xemiss + xw + xq + xg;
721   Float_t ytot = yemiss + yw + yq + yg;
722
723   SetXPointOnCathode(xtot);
724   SetYPointOnCathode(ytot);
725
726
727   Float_t distanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2)
728                                     +TMath::Power(fPhotonLimitY,2)); 
729
730   return distanceFromEntrance;
731
732 }
733 //__________________________________________________________________________________________________
734 void AliRICHRecon::FindPhiPoint()
735 {
736   //find phi of generated point 
737
738   Float_t xtoentr = GetEntranceX();
739   Float_t ytoentr = GetEntranceY();
740
741   Float_t trackTheta = GetTrackTheta();
742   Float_t trackPhi = GetTrackPhi();
743
744   Float_t emissionPoint = GetEmissionPoint();
745
746   Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi);
747   Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi);
748   Float_t phi = atan2(argY,argX);
749
750   SetPhiPoint(phi);
751
752 }
753 //__________________________________________________________________________________________________
754 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
755 {
756   // cerenkov angle from n and beta
757
758 // Compute the cerenkov angle
759
760   Float_t thetacer;
761
762   if((n*beta)<1.) {
763     thetacer = 999.;
764     //    cout << " warning in Cerenkoangle !!!!!! " << endl;
765     return thetacer;
766   }
767
768   thetacer = acos (1./(n*beta));
769   return thetacer;
770 }
771 //__________________________________________________________________________________________________
772 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
773 {
774   // Snell law
775
776 // Compute the Snell angle
777
778   Float_t sinrefractangle;
779   Float_t refractangle;
780
781   sinrefractangle = (n1/n2)*sin(theta1);
782
783   if(sinrefractangle>1.) {
784     //    cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
785     refractangle = 999.;
786     return refractangle;
787   }
788   
789   refractangle = asin(sinrefractangle);  
790   return refractangle;
791 }
792 //__________________________________________________________________________________________________
793 void AliRICHRecon::HoughResponse()
794 {
795   //Hough response
796
797 // Implement Hough response pat. rec. method
798
799   Float_t *hCSspace;
800
801   int           bin=0;
802   int           bin1=0;
803   int           bin2=0;
804   int           i, j, k, nCorrBand;
805   float         hcs[750],hcsw[750];
806   float         angle, weight;
807   float         lowerlimit,upperlimit;
808
809   float         etaPeak[100];
810
811   int           nBin;
812
813   float etaPeakPos  = -1;
814
815   Int_t   etaPeakCount = -1;
816   
817   Float_t thetaCerenkov = 0.;
818     
819   nBin = (int)(0.5+fThetaMax/(fDTheta));
820   nCorrBand = (int)(0.5+ fWindowWidth/(2 * fDTheta)); 
821
822   memset ((void *)hcs, 0, fThetaBin*sizeof(float));
823   memset ((void *)hcsw, 0, fThetaBin*sizeof(float));
824
825   Int_t nPhotons = GetPhotonsNumber();
826
827   Int_t weightFlag = 0;
828
829   for (k=0; k< nPhotons; k++) {
830
831     SetPhotonIndex(k);
832
833     angle = GetPhotonEta();
834
835     if(angle == -999.) continue;
836
837     if (angle>=fThetaMin && angle<= fThetaMax) 
838
839       {
840
841         bin = (int)(0.5+angle/(fDTheta));
842
843         bin1= bin-nCorrBand;
844         bin2= bin+nCorrBand;
845
846         // calculate weights
847
848         if(fIsWEIGHT)
849           {
850             lowerlimit = ((Float_t)bin1)*fDTheta + 0.5*fDTheta;
851             SetThetaOfRing(lowerlimit);
852             FindAreaAndPortionOfRing();
853             Float_t area1 = GetAreaOfRing();
854             
855             upperlimit = ((Float_t)bin2)*fDTheta + 0.5*fDTheta;
856             SetThetaOfRing(upperlimit);
857             FindAreaAndPortionOfRing();
858             Float_t area2 = GetAreaOfRing();
859             
860             //      cout << "lowerlimit" << lowerlimit << "upperlimit " << upperlimit << endl;
861             Float_t diffarea = area2 - area1;
862
863             if(diffarea>0)
864               {
865                 weight = 1./(area2-area1);
866               }
867             else
868               {
869                 weightFlag = 1;
870                 weight = 1.;
871               }
872
873             //      cout <<" low "<< lowerlimit << " up " << upperlimit << 
874             //        " area1 " << area1 << " area2 " << area2 << " weight " << weight << endl;
875             
876           }
877         else
878           {
879             weight = 1.;
880           }
881
882         SetPhotonWeight(weight);
883         
884         //      cout << "weight..." << weight << endl;
885
886
887         if (bin1<0)    bin1=0;
888         if (bin2>nBin) bin2=nBin;
889       
890         for (j=bin1; j<bin2; j++) 
891           {
892             hcs[j] += 1; 
893             hcsw[j] += weight;
894           }
895       }
896   }
897   
898
899   if(weightFlag == 0) 
900     {
901       hCSspace = hcsw;
902     }
903   else
904     {
905       hCSspace = hcs;
906       //      cout << " probems with weight...normal procedure adopted " << endl;
907     }
908
909   HoughFiltering(hCSspace);
910
911   for (bin=0; bin <nBin; bin++) {
912     angle = (bin+0.5) * (fDTheta);
913     if (hCSspace[bin] && hCSspace[bin] > etaPeakPos) {
914       etaPeakCount = 0;
915       etaPeakPos = hCSspace[bin];
916       etaPeak[0]=angle;
917     }
918     else { 
919       if (hCSspace[bin] == etaPeakPos) {
920         etaPeak[++etaPeakCount] = angle;
921       }
922     }
923   } 
924
925   for (i=0; i<etaPeakCount+1; i++) {
926     thetaCerenkov += etaPeak[i];
927   }
928   if (etaPeakCount>=0) {
929     thetaCerenkov /= etaPeakCount+1;
930     fThetaPeakPos = etaPeakPos;
931   }
932
933   SetThetaCerenkov(thetaCerenkov);
934 }
935 //__________________________________________________________________________________________________
936 void AliRICHRecon::HoughFiltering(float hcs[])
937 {
938   // filter for Hough
939
940 // hough filtering
941
942    float hcsFilt[750];
943    float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
944    int nx, i, nxDx;
945    int sizeHCS;
946    int nBin;
947
948    nBin =  (int)(1+fThetaMax/fDTheta); 
949    sizeHCS = fThetaBin*sizeof(float);
950
951    memset ((void *)hcsFilt, 0, sizeHCS); 
952
953    for (nx = 0; nx < nBin; nx++) {
954       for (i = 0; i < 5; i++)   {
955         nxDx = nx + (i-2);
956         if (nxDx> -1 && nxDx<nBin)
957              hcsFilt[nx] +=  hcs[nxDx] * k[i];
958       }      
959    }
960      
961    for (nx = 0; nx < nBin; nx++) {
962      hcs[nx] = hcsFilt[nx];
963    }
964 }
965 //__________________________________________________________________________________________________
966 void AliRICHRecon::FindThetaCerenkov()
967 {
968   // manage with weight for photons
969
970   Float_t wei = 0.;
971   Float_t weightThetaCerenkov = 0.;
972
973   Int_t nPhotons = GetPhotonsNumber();
974   for(Int_t i=0;i<nPhotons;i++)
975     {
976       SetPhotonIndex(i);
977
978       if(GetPhotonFlag() == 2)
979         {
980           Float_t photonEta = GetPhotonEta();
981           Float_t photonWeight = GetPhotonWeight();
982           weightThetaCerenkov += photonEta*photonWeight;
983           wei += photonWeight;
984         }
985     }
986
987   if(wei != 0.) 
988     {
989       weightThetaCerenkov /= wei;
990     }
991   else
992     {
993       weightThetaCerenkov = 0.;
994     }
995   
996   SetThetaCerenkov(weightThetaCerenkov);
997
998   AliDebug(1,Form(" thetac weighted -> %f",weightThetaCerenkov));
999 }
1000 //__________________________________________________________________________________________________
1001 void AliRICHRecon::FlagPhotons()
1002 {
1003   // flag photons
1004
1005   Int_t nPhotonHough = 0;
1006
1007   Float_t thetaCerenkov = GetThetaCerenkov();
1008   AliDebug(1,Form(" fThetaCerenkov %f ",thetaCerenkov));
1009
1010   Float_t thetaDist= thetaCerenkov - fThetaMin;
1011   Int_t steps = (Int_t)(thetaDist / fDTheta);
1012
1013   Float_t tmin = fThetaMin + (Float_t)(steps - 1)*fDTheta;
1014   Float_t tmax = fThetaMin + (Float_t)(steps)*fDTheta;
1015   Float_t tavg = 0.5*(tmin+tmax);
1016
1017   tmin = tavg - 0.5*fWindowWidth;
1018   tmax = tavg + 0.5*fWindowWidth;
1019
1020   //  Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
1021
1022   Int_t nPhotons = GetPhotonsNumber();
1023
1024   //  for(Int_t i=0;i<candidatePhotonsNumber;i++)
1025
1026   for(Int_t i=0;i<nPhotons;i++)
1027     {
1028       SetPhotonIndex(i);
1029
1030       Float_t photonEta = GetPhotonEta();
1031
1032       if(photonEta == -999.) continue;
1033
1034       if(photonEta >= tmin && photonEta <= tmax)
1035         {
1036           SetPhotonFlag(2);
1037           nPhotonHough++;
1038         }
1039     }
1040   SetHoughPhotons(nPhotonHough);
1041 }
1042