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