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