RICH/AliRICHPatRec.cxx

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