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