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