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