5ceada666dca8df06e922c02da8a012c1338e32b
[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.060;
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 kBad;//no clusters at all for a given track
67   Bool_t kPatRec = kFALSE;  
68
69   Int_t candidatePhotons = 0;
70
71   SetThetaCerenkov(999.);
72   SetHoughPhotons(0);
73   SetHoughPhotonsNorm(0);
74
75   for (Int_t j=0; j < fpClusters->GetEntries(); j++){//clusters loop
76     SetPhotonIndex(j);
77     SetPhotonFlag(0);
78     SetPhotonEta(-999.);
79     SetPhotonWeight(0.);
80     if (j == GetMipIndex()) continue; // do not consider MIP cluster as a candidate photon
81     Float_t xtoentr = ((AliRICHcluster*)fpClusters->UncheckedAt(j))->X() - GetShiftX();
82     Float_t ytoentr = ((AliRICHcluster*)fpClusters->UncheckedAt(j))->Y() - GetShiftY();
83     SetEntranceX(xtoentr);
84     SetEntranceY(ytoentr);
85     FindPhiPoint();
86 //      Int_t photonStatus = PhotonInBand();
87 //      if(photonStatus == 0) continue;
88     SetPhotonFlag(1);
89     FindThetaPhotonCerenkov();
90     Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
91     AliDebug(1,Form("THETA CERENKOV ---> %i",thetaPhotonCerenkov));
92     SetPhotonEta(thetaPhotonCerenkov);
93     candidatePhotons++;
94   }//clusters loop
95
96   if(candidatePhotons >= 1) kPatRec = kTRUE;
97
98   if(!kPatRec) return kBad;
99
100   SetPhotonsNumber(fpClusters->GetEntries());
101
102   HoughResponse();
103   
104   fNrings++;
105
106   FlagPhotons();
107   Int_t nPhotonHough = GetHoughPhotons();
108  
109   if(nPhotonHough < 1) 
110     {
111       SetThetaCerenkov(999.);
112       SetHoughPhotonsNorm(0.);
113       return kBad;
114     }
115
116   if(fIsWEIGHT) FindWeightThetaCerenkov();
117
118   Float_t thetaCerenkov = GetThetaCerenkov();
119
120   AliDebug(1,Form("Number of clusters accepted --->  %i",nPhotonHough));
121   
122   SetThetaOfRing(thetaCerenkov);
123 //  FindAreaAndPortionOfRing();
124
125 //  Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
126 //  SetHoughPhotonsNorm(nPhotonHoughNorm);
127
128   // Calculate the area where the photon are accepted...
129 /*
130   Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth; 
131   SetThetaOfRing(thetaInternal);
132   FindAreaAndPortionOfRing();
133   Float_t internalArea = GetAreaOfRing();
134
135   Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth; 
136   SetThetaOfRing(thetaExternal);
137   FindAreaAndPortionOfRing();
138   Float_t externalArea = GetAreaOfRing();
139
140   Float_t houghArea = externalArea - internalArea;
141
142   SetHoughArea(houghArea);
143 */
144   return GetThetaCerenkov();
145
146 }//ThetaCerenkov()
147 //__________________________________________________________________________________________________
148 void AliRICHRecon::FindEmissionPoint()
149 {
150   //estimate the emission point in radiator
151
152 // Find emission point
153
154   Float_t absorbtionLenght=7.83*fRadiatorWidth; //absorption length in the freon (cm)
155   // 7.83 = -1/ln(T0) where 
156   // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
157   Float_t photonLenght, photonLenghtMin, photonLenghtMax;
158
159   photonLenght=exp(-fRadiatorWidth/(absorbtionLenght*cos(fCerenkovAnglePad)));
160   photonLenghtMin=fRadiatorWidth*photonLenght/(1.-photonLenght);
161   photonLenghtMax=absorbtionLenght*cos(fCerenkovAnglePad);
162   Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax;
163
164   SetEmissionPoint(emissionPoint);
165   SetEmissionPoint(fRadiatorWidth/2); // tune the emission point
166 }
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(6.85);
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(6.85);
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
735 void AliRICHRecon::FindPhiPoint()
736 {
737   //find phi of generated point 
738
739   Float_t xtoentr = GetEntranceX();
740   Float_t ytoentr = GetEntranceY();
741
742   Float_t trackTheta = GetTrackTheta();
743   Float_t trackPhi = GetTrackPhi();
744
745   Float_t emissionPoint = GetEmissionPoint();
746
747   Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi);
748   Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi);
749   Float_t phi = atan2(argY,argX); 
750
751   SetPhiPoint(phi);
752
753 }
754
755 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
756 {
757   // cerenkov angle from n and beta
758
759 // Compute the cerenkov angle
760
761   Float_t thetacer;
762
763   if((n*beta)<1.) {
764     thetacer = 999.;
765     //    cout << " warning in Cerenkoangle !!!!!! " << endl;
766     return thetacer;
767   }
768
769   thetacer = acos (1./(n*beta));
770   return thetacer;
771 }
772
773 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
774 {
775   // Snell law
776
777 // Compute the Snell angle
778
779   Float_t sinrefractangle;
780   Float_t refractangle;
781
782   sinrefractangle = (n1/n2)*sin(theta1);
783
784   if(sinrefractangle>1.) {
785     //    cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
786     refractangle = 999.;
787     return refractangle;
788   }
789   
790   refractangle = asin(sinrefractangle);  
791   return refractangle;
792 }
793
794
795 void AliRICHRecon::HoughResponse()
796 {
797   //Hough response
798
799 // Implement Hough response pat. rec. method
800
801   Float_t *hCSspace;
802
803   int           bin=0;
804   int           bin1=0;
805   int           bin2=0;
806   int           i, j, k, nCorrBand;
807   float         hcs[750],hcsw[750];
808   float         angle, weight;
809   float         lowerlimit,upperlimit;
810
811   float         etaPeak[100];
812
813   int           nBin;
814
815   float etaPeakPos  = -1;
816
817   Int_t   etaPeakCount = -1;
818   
819   Float_t thetaCerenkov = 0.;
820     
821   nBin = (int)(0.5+fThetaMax/(fDTheta));
822   nCorrBand = (int)(0.5+ fWindowWidth/(2 * fDTheta)); 
823
824   memset ((void *)hcs, 0, fThetaBin*sizeof(float));
825   memset ((void *)hcsw, 0, fThetaBin*sizeof(float));
826
827   Int_t nPhotons = GetPhotonsNumber();
828
829   Int_t weightFlag = 0;
830
831   for (k=0; k< nPhotons; k++) {
832
833     SetPhotonIndex(k);
834
835     angle = GetPhotonEta();
836
837     if(angle == -999.) continue;
838
839     if (angle>=fThetaMin && angle<= fThetaMax) 
840
841       {
842
843         bin = (int)(0.5+angle/(fDTheta));
844
845         bin1= bin-nCorrBand;
846         bin2= bin+nCorrBand;
847
848         // calculate weights
849
850         if(fIsWEIGHT)
851           {
852             lowerlimit = ((Float_t)bin1)*fDTheta + 0.5*fDTheta;
853             SetThetaOfRing(lowerlimit);
854             FindAreaAndPortionOfRing();
855             Float_t area1 = GetAreaOfRing();
856             
857             upperlimit = ((Float_t)bin2)*fDTheta + 0.5*fDTheta;
858             SetThetaOfRing(upperlimit);
859             FindAreaAndPortionOfRing();
860             Float_t area2 = GetAreaOfRing();
861             
862             //      cout << "lowerlimit" << lowerlimit << "upperlimit " << upperlimit << endl;
863             Float_t diffarea = area2 - area1;
864
865             if(diffarea>0)
866               {
867                 weight = 1./(area2-area1);
868               }
869             else
870               {
871                 weightFlag = 1;
872                 weight = 1.;
873               }
874
875             //      cout <<" low "<< lowerlimit << " up " << upperlimit << 
876             //        " area1 " << area1 << " area2 " << area2 << " weight " << weight << endl;
877             
878           }
879         else
880           {
881             weight = 1.;
882           }
883
884         SetPhotonWeight(weight);
885         
886         //      cout << "weight..." << weight << endl;
887
888
889         if (bin1<0)    bin1=0;
890         if (bin2>nBin) bin2=nBin;
891       
892         for (j=bin1; j<bin2; j++) 
893           {
894             hcs[j] += 1; 
895             hcsw[j] += weight;
896           }
897       }
898   }
899   
900
901   if(weightFlag == 0) 
902     {
903       hCSspace = hcsw;
904     }
905   else
906     {
907       hCSspace = hcs;
908       //      cout << " probems with weight...normal procedure adopted " << endl;
909     }
910
911   HoughFiltering(hCSspace);
912
913   for (bin=0; bin <nBin; bin++) {
914     angle = (bin+0.5) * (fDTheta);
915     if (hCSspace[bin] && hCSspace[bin] > etaPeakPos) {
916       etaPeakCount = 0;
917       etaPeakPos = hCSspace[bin];
918       etaPeak[0]=angle;
919     }
920     else { 
921       if (hCSspace[bin] == etaPeakPos) {
922         etaPeak[++etaPeakCount] = angle;
923       }
924     }
925   } 
926
927   for (i=0; i<etaPeakCount+1; i++) {
928     thetaCerenkov += etaPeak[i];
929   }
930   if (etaPeakCount>=0) {
931     thetaCerenkov /= etaPeakCount+1;
932     fThetaPeakPos = etaPeakPos;
933   }
934
935   SetThetaCerenkov(thetaCerenkov);
936 }
937
938
939 void AliRICHRecon::HoughFiltering(float hcs[])
940 {
941   // filter for Hough
942
943 // hough filtering
944
945    float hcsFilt[750];
946    float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
947    int nx, i, nxDx;
948    int sizeHCS;
949    int nBin;
950
951    nBin =  (int)(1+fThetaMax/fDTheta); 
952    sizeHCS = fThetaBin*sizeof(float);
953
954    memset ((void *)hcsFilt, 0, sizeHCS); 
955
956    for (nx = 0; nx < nBin; nx++) {
957       for (i = 0; i < 5; i++)   {
958         nxDx = nx + (i-2);
959         if (nxDx> -1 && nxDx<nBin)
960              hcsFilt[nx] +=  hcs[nxDx] * k[i];
961       }      
962    }
963      
964    for (nx = 0; nx < nBin; nx++) {
965      hcs[nx] = hcsFilt[nx];
966    }
967 }
968
969 void AliRICHRecon::FindWeightThetaCerenkov()
970 {
971   // manage with weight for photons
972
973   Float_t wei = 0.;
974   Float_t weightThetaCerenkov = 0.;
975
976   Int_t nPhotons = GetPhotonsNumber();
977   for(Int_t i=0;i<nPhotons;i++)
978     {
979       SetPhotonIndex(i);
980
981       if(GetPhotonFlag() == 2)
982         {
983           Float_t photonEta = GetPhotonEta();
984           Float_t photonWeight = GetPhotonWeight();
985           weightThetaCerenkov += photonEta*photonWeight;
986           wei += photonWeight;
987         }
988     }
989
990   if(wei != 0.) 
991     {
992       weightThetaCerenkov /= wei;
993     }
994   else
995     {
996       weightThetaCerenkov = 0.;
997     }
998   
999   SetThetaCerenkov(weightThetaCerenkov);
1000
1001   cout << " thetac weighted -> " << weightThetaCerenkov << endl;
1002 }
1003
1004
1005 void AliRICHRecon::FlagPhotons()
1006 {
1007   // flag photons
1008
1009   Int_t nPhotonHough = 0;
1010
1011   Float_t thetaCerenkov = GetThetaCerenkov();
1012   AliDebug(1,Form(" fThetaCerenkov ",thetaCerenkov));
1013
1014   Float_t thetaDist= thetaCerenkov - fThetaMin;
1015   Int_t steps = (Int_t)(thetaDist / fDTheta);
1016
1017   Float_t tmin = fThetaMin + (Float_t)(steps - 1)*fDTheta;
1018   Float_t tmax = fThetaMin + (Float_t)(steps)*fDTheta;
1019   Float_t tavg = 0.5*(tmin+tmax);
1020
1021   tmin = tavg - 0.5*fWindowWidth;
1022   tmax = tavg + 0.5*fWindowWidth;
1023
1024   //  Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
1025
1026   Int_t nPhotons = GetPhotonsNumber();
1027
1028   //  for(Int_t i=0;i<candidatePhotonsNumber;i++)
1029
1030   for(Int_t i=0;i<nPhotons;i++)
1031     {
1032       SetPhotonIndex(i);
1033
1034       Float_t photonEta = GetPhotonEta();
1035
1036       if(photonEta == -999.) continue;
1037
1038       if(photonEta >= tmin && photonEta <= tmax)
1039         {
1040           SetPhotonFlag(2);
1041           nPhotonHough++;
1042         }
1043     }
1044   SetHoughPhotons(nPhotonHough);
1045 }
1046