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