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