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