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