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