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