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