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