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