RICH/AliRICHRecon.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 //                                                                      //
  18 // AliRICHRecon                                                         //
  19 //                                                                      //
  20 // RICH class to perfom pattern recognition based on Hough transfrom    //
  21 //                                                                      //
  22 //////////////////////////////////////////////////////////////////////////
  23
  24 #include <Riostream.h>
  25 #include <TMath.h>
  26 #include <TRotation.h>
  27 #include <TVector3.h>
  28
  29 #include "AliRICH.h"
  30 #include "AliRICHParam.h"
  31 #include "AliRICHRecon.h"
  32 #include "AliRICHHelix.h"
  33 #include <AliLog.h>
  34
  35 #define NPointsOfRing 201
  36
  37 //__________________________________________________________________________________________________
  38 AliRICHRecon::AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId)
  39              :TTask("RichRec","RichPat")
  40 {
  41 // main ctor
  42   SetFreonScaleFactor(1);
  43   fIsWEIGHT = kFALSE;
  44   fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75;
  45   fDTheta       = 0.001;   fWindowWidth  = 0.060;
  46   fRadiatorWidth = AliRICHParam::Zfreon();
  47   fQuartzWidth   = AliRICHParam::Zwin();
  48   fGapWidth      = AliRICHParam::Freon2Pc() - fRadiatorWidth - fQuartzWidth;
  49   fXmin         = -AliRICHParam::PcSizeX()/2.;
  50   fXmax         =  AliRICHParam::PcSizeX()/2.;
  51   fYmin         = -AliRICHParam::PcSizeY()/2.;
  52   fYmax         =  AliRICHParam::PcSizeY()/2.;
  53   SetTrackTheta(pHelix->Ploc().Theta());
  54   SetTrackPhi(pHelix->Ploc().Phi());
  55   SetMipIndex(iMipId);
  56   SetShiftX(pHelix->PosRad().X());
  57   SetShiftY(pHelix->PosRad().Y());
  58   fpClusters = pClusters;
  59 }
  60 //__________________________________________________________________________________________________
  61 Double_t AliRICHRecon::ThetaCerenkov()
  62 {
  63 // Pattern recognition method based on Hough transform
  64 // Return theta Cerenkov for a given track and list of clusters which are set in ctor
  65
  66   if(fpClusters->GetEntries()==0) return kBad;//no clusters at all for a given track
  67   Bool_t kPatRec = kFALSE;
  68
  69   AliDebug(1,Form("---Track Parameters--- Theta: %f , Phi: %f ",GetTrackTheta()*TMath::RadToDeg(),GetTrackPhi()*TMath::RadToDeg()));
  70
  71   Int_t candidatePhotons = 0;
  72
  73   SetThetaCerenkov(999.);
  74   SetHoughPhotons(0);
  75   SetHoughPhotonsNorm(0);
  76
  77   for (Int_t j=0; j < fpClusters->GetEntries(); j++){//clusters loop
  78     SetPhotonIndex(j);
  79     SetPhotonFlag(0);
  80     SetPhotonEta(-999.);
  81     SetPhotonWeight(0.);
  82     if (j == GetMipIndex()) continue; // do not consider MIP cluster as a candidate photon
  83     Float_t xtoentr = ((AliRICHcluster*)fpClusters->UncheckedAt(j))->X() - GetShiftX();
  84     Float_t ytoentr = ((AliRICHcluster*)fpClusters->UncheckedAt(j))->Y() - GetShiftY();
  85     SetEntranceX(xtoentr);
  86     SetEntranceY(ytoentr);
  87     FindPhiPoint();
  88 //      Int_t photonStatus = PhotonInBand();
  89 //      if(photonStatus == 0) continue;
  90     SetPhotonFlag(1);
  91     FindThetaPhotonCerenkov();
  92     Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
  93     AliDebug(1,Form("THETA CERENKOV ---> %f",thetaPhotonCerenkov));
  94     SetPhotonEta(thetaPhotonCerenkov);
  95     candidatePhotons++;
  96   }//clusters loop
  97
  98   if(candidatePhotons >= 1) kPatRec = kTRUE;
  99
 100   if(!kPatRec) return kBad;
 101
 102   SetPhotonsNumber(fpClusters->GetEntries());
 103
 104   HoughResponse();
 105
 106   fNrings++;
 107
 108   FlagPhotons();
 109   Int_t nPhotonHough = GetHoughPhotons();
 110
 111   if(nPhotonHough < 1)
 112     {
 113       SetThetaCerenkov(999.);
 114       SetHoughPhotonsNorm(0.);
 115       return kBad;
 116     }
 117
 118   if(fIsWEIGHT) FindWeightThetaCerenkov();
 119
 120   Float_t thetaCerenkov = GetThetaCerenkov();
 121
 122   AliDebug(1,Form("Number of clusters accepted --->  %i",nPhotonHough));
 123
 124   SetThetaOfRing(thetaCerenkov);
 125 //  FindAreaAndPortionOfRing();
 126
 127 //  Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
 128 //  SetHoughPhotonsNorm(nPhotonHoughNorm);
 129
 130   // Calculate the area where the photon are accepted...
 131 /*
 132   Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth;
 133   SetThetaOfRing(thetaInternal);
 134   FindAreaAndPortionOfRing();
 135   Float_t internalArea = GetAreaOfRing();
 136
 137   Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth;
 138   SetThetaOfRing(thetaExternal);
 139   FindAreaAndPortionOfRing();
 140   Float_t externalArea = GetAreaOfRing();
 141
 142   Float_t houghArea = externalArea - internalArea;
 143
 144   SetHoughArea(houghArea);
 145 */
 146   return GetThetaCerenkov();
 147
 148 }//ThetaCerenkov()
 149 //__________________________________________________________________________________________________
 150 void AliRICHRecon::FindEmissionPoint()
 151 {
 152   //estimate the emission point in radiator
 153
 154 // Find emission point
 155
 156   Float_t absorbtionLenght=7.83*fRadiatorWidth; //absorption length in the freon (cm)
 157   // 7.83 = -1/ln(T0) where
 158   // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
 159   Float_t photonLenght, photonLenghtMin, photonLenghtMax;
 160
 161   photonLenght=exp(-fRadiatorWidth/(absorbtionLenght*cos(fCerenkovAnglePad)));
 162   photonLenghtMin=fRadiatorWidth*photonLenght/(1.-photonLenght);
 163   photonLenghtMax=absorbtionLenght*cos(fCerenkovAnglePad);
 164   Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax;
 165
 166   SetEmissionPoint(emissionPoint);
 167   SetEmissionPoint(fRadiatorWidth/2); // tune the emission point
 168 }
 169
 170
 171 Int_t AliRICHRecon::PhotonInBand()
 172 {
 173   //search band fro photon candidates
 174
 175   //  Float_t massOfParticle;
 176   Float_t nfreon;
 177
 178   Float_t thetacer;
 179
 180   Float_t xtoentr = GetEntranceX();
 181   Float_t ytoentr = GetEntranceY();
 182
 183   Float_t innerRadius;
 184   Float_t outerRadius;
 185
 186   Float_t phpad = GetPhiPoint();
 187
 188
 189   // inner radius //
 190   SetPhotonEnergy(5.6);
 191   SetEmissionPoint(fRadiatorWidth -0.0001);
 192   SetFreonRefractiveIndex();
 193
 194   nfreon = GetFreonRefractiveIndex();
 195   thetacer = 0.;
 196
 197   AliDebug(1,Form("thetacer in photoninband min %f",thetacer));
 198
 199   FindThetaAtQuartz(thetacer);
 200
 201   if(thetacer == 999. || GetThetaAtQuartz() == 999.)
 202     {
 203       innerRadius = -999.;
 204       SetXInnerRing(-999.);
 205       SetYInnerRing(-999.);
 206       SetRadiusInnerRing(-999.);
 207     }
 208   else
 209     {
 210       SetThetaPhotonInDRS(GetThetaAtQuartz());
 211       SetPhiPhotonInDRS(phpad);
 212
 213       innerRadius = FromEmissionToCathode();
 214        if(innerRadius == 999.) innerRadius = -999.;
 215
 216       SetXInnerRing(GetXPointOnCathode());
 217       SetYInnerRing(GetYPointOnCathode());
 218       SetRadiusInnerRing(innerRadius);
 219     }
 220
 221   // outer radius //
 222   SetPhotonEnergy(7.7);
 223   SetEmissionPoint(0.);
 224 //  SetMassHypotesis(0.139567);
 225   SetFreonRefractiveIndex();
 226
 227   nfreon = GetFreonRefractiveIndex();
 228
 229   thetacer = Cerenkovangle(nfreon,1);
 230
 231   //  thetacer = 0.75;
 232
 233   AliDebug(1,Form("thetacer in photoninband max %f",thetacer));
 234
 235   FindThetaAtQuartz(thetacer);
 236
 237   if(thetacer == 999. || GetThetaAtQuartz() == 999.)
 238     {
 239       outerRadius = 999.;
 240       SetXOuterRing(999.);
 241       SetYOuterRing(999.);
 242       SetRadiusOuterRing(999.);
 243     }
 244   else
 245     {
 246       SetThetaPhotonInDRS(GetThetaAtQuartz());
 247       SetPhiPhotonInDRS(phpad);
 248
 249       outerRadius = FromEmissionToCathode();
 250 //      cout << " outerRadius " << outerRadius << endl;
 251       SetXOuterRing(GetXPointOnCathode());
 252       SetYOuterRing(GetYPointOnCathode());
 253       SetRadiusOuterRing(outerRadius);
 254     }
 255
 256   Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2));
 257
 258   AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius));
 259
 260   if(padradius>=innerRadius && padradius<=outerRadius) return 1;
 261   return 0;
 262 }
 263
 264 void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
 265 {
 266   //find the theta at the quartz plate
 267
 268   if(thetaCerenkov == 999.)
 269     {
 270       SetThetaAtQuartz(999.);
 271       return;
 272     }
 273
 274   Float_t thetaAtQuartz = 999.;
 275
 276   Float_t trackTheta = GetTrackTheta();
 277
 278   if(trackTheta == 0) {
 279
 280     thetaAtQuartz = thetaCerenkov;
 281     SetThetaAtQuartz(thetaAtQuartz);
 282     return;
 283   }
 284
 285   Float_t trackPhi   = GetTrackPhi();
 286   Float_t phiPoint = GetPhiPoint();
 287
 288   Double_t den = TMath::Sin((Double_t)trackTheta)
 289     *TMath::Cos((Double_t)trackPhi)
 290     *TMath::Cos((Double_t)phiPoint) +
 291     TMath::Sin((Double_t)trackTheta)
 292     *TMath::Sin((Double_t)trackPhi)
 293     *TMath::Sin((Double_t)phiPoint);
 294   Double_t b = TMath::Cos((Double_t)trackTheta)/den;
 295   Double_t c = -TMath::Cos((Double_t)thetaCerenkov)/den;
 296
 297   Double_t underSqrt = 1 + b*b - c*c;
 298
 299   if(underSqrt < 0) {
 300     SetThetaAtQuartz(999.);
 301     return;
 302   }
 303
 304   Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c);
 305   Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c);
 306
 307   Double_t thetaSol1 = 2*TMath::ATan(sol1);
 308   Double_t thetaSol2 = 2*TMath::ATan(sol2);
 309
 310   if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1;
 311   if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2;
 312
 313 //  AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz));
 314
 315   SetThetaAtQuartz(thetaAtQuartz);
 316 }
 317
 318 void AliRICHRecon::FindThetaPhotonCerenkov()
 319 {
 320   //find theta cerenkov of ring
 321
 322   Float_t thetaCerMin = 0.;
 323   Float_t thetaCerMax = 0.75;
 324   Float_t thetaCerMean;
 325
 326   Float_t radiusMin, radiusMax, radiusMean;
 327   Int_t nIteration = 0;
 328
 329   const Float_t kTollerance = 0.05;
 330
 331
 332   Float_t phiPoint = GetPhiPoint();
 333
 334   SetPhotonEnergy(6.85);
 335   SetEmissionPoint(fRadiatorWidth/2);
 336
 337   Float_t xPoint = GetEntranceX();
 338   Float_t yPoint = GetEntranceY();
 339   Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint);
 340
 341   AliDebug(1,Form(" DistPoint %f ",distPoint));
 342
 343   // Star minimization...
 344
 345   // First value...
 346
 347   FindThetaAtQuartz(thetaCerMin);
 348
 349   if(GetThetaAtQuartz() == 999.)
 350     {
 351       radiusMin = -999.;
 352     }
 353   else
 354     {
 355       SetThetaPhotonInDRS(GetThetaAtQuartz());
 356       SetPhiPhotonInDRS(phiPoint);
 357
 358       radiusMin = FromEmissionToCathode();
 359     }
 360
 361   // Second value...
 362
 363   FindThetaAtQuartz(thetaCerMax);
 364   if(GetThetaAtQuartz() == 999.)
 365     {
 366       radiusMax = 999.;
 367     }
 368   else
 369     {
 370       SetThetaPhotonInDRS(GetThetaAtQuartz());
 371       SetPhiPhotonInDRS(phiPoint);
 372
 373       radiusMax = FromEmissionToCathode();
 374     }
 375   // Mean value...
 376
 377   thetaCerMean = (thetaCerMax + thetaCerMin)/2;
 378
 379   FindThetaAtQuartz(thetaCerMean);
 380   if(GetThetaAtQuartz() == 999.)
 381     {
 382       radiusMean = 999.;
 383     }
 384   else
 385     {
 386       SetThetaPhotonInDRS(GetThetaAtQuartz());
 387       SetPhiPhotonInDRS(phiPoint);
 388
 389       radiusMean = FromEmissionToCathode();
 390     }
 391
 392   AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax));
 393
 394   while (TMath::Abs(radiusMean-distPoint) > kTollerance)
 395     {
 396
 397       if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean;
 398       if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) {
 399
 400         thetaCerMin = thetaCerMean;
 401
 402         FindThetaAtQuartz(thetaCerMin);
 403         SetThetaPhotonInDRS(GetThetaAtQuartz());
 404         SetPhiPhotonInDRS(phiPoint);
 405
 406         radiusMin =FromEmissionToCathode();
 407       }
 408
 409       thetaCerMean = (thetaCerMax + thetaCerMin)/2;
 410
 411       FindThetaAtQuartz(thetaCerMean);
 412       SetThetaPhotonInDRS(GetThetaAtQuartz());
 413       SetPhiPhotonInDRS(phiPoint);
 414
 415       radiusMean = FromEmissionToCathode();
 416
 417       nIteration++;
 418       if(nIteration>=50) {
 419         AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration));
 420         SetThetaPhotonCerenkov(999.);
 421         return;
 422       }
 423     }
 424
 425   AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean));
 426   SetThetaPhotonCerenkov(thetaCerMean);
 427
 428 }
 429
 430 void AliRICHRecon::FindAreaAndPortionOfRing()
 431 {
 432   //find fraction of the ring accepted by the RICH
 433
 434   Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing];
 435
 436   //  Float_t xtoentr = GetEntranceX();
 437   //  Float_t ytoentr = GetEntranceY();
 438   Float_t shiftX = GetShiftX();
 439   Float_t shiftY = GetShiftY();
 440
 441   Float_t xemiss = GetXCoordOfEmission();
 442   Float_t yemiss = GetYCoordOfEmission();
 443
 444   Float_t x0 = xemiss + shiftX;
 445   Float_t y0 = yemiss + shiftY;
 446
 447
 448   SetPhotonEnergy(6.85);
 449   SetFreonRefractiveIndex();
 450
 451   SetEmissionPoint(fRadiatorWidth/2.);
 452
 453   Float_t theta = GetThetaOfRing();
 454
 455   Int_t nPoints = 0;
 456   Int_t nPsiAccepted = 0;
 457   Int_t nPsiTotal = 0;
 458
 459   for(Int_t i=0;i<NPointsOfRing-1;i++)
 460     {
 461
 462       Float_t psi = 2*TMath::Pi()*i/NPointsOfRing;
 463
 464       SetThetaPhotonInTRS(theta);
 465       SetPhiPhotonInTRS(psi);
 466       FindPhotonAnglesInDRS();
 467
 468       Float_t radius = FromEmissionToCathode();
 469       if (radius == 999.) continue;
 470
 471       nPsiTotal++;
 472
 473       Float_t xPointRing = GetXPointOnCathode() + shiftX;
 474       Float_t yPointRing = GetYPointOnCathode() + shiftY;
 475
 476       SetDetectorWhereX(xPointRing);
 477       SetDetectorWhereY(yPointRing);
 478
 479       Int_t zone = CheckDetectorAcceptance();
 480
 481
 482       if (zone != 0)
 483         {
 484           FindIntersectionWithDetector();
 485           xPoint[nPoints] = GetIntersectionX();
 486           yPoint[nPoints] = GetIntersectionY();
 487         }
 488       else
 489         {
 490           xPoint[nPoints] = xPointRing;
 491           yPoint[nPoints] = yPointRing;
 492           nPsiAccepted++;
 493         }
 494
 495       nPoints++;
 496
 497     }
 498
 499   xPoint[nPoints] = xPoint[0];
 500   yPoint[nPoints] = yPoint[0];
 501
 502   // find area...
 503
 504   Float_t area = 0;
 505
 506   for (Int_t i = 0; i < nPoints; i++)
 507     {
 508       area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0));
 509     }
 510
 511   area *= 0.5;
 512
 513   Float_t portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal));
 514
 515
 516   SetAreaOfRing(area);
 517   SetPortionOfRing(portionOfRing);
 518 }
 519
 520 void AliRICHRecon::FindIntersectionWithDetector()
 521 {
 522   // find ring intersection with CsI edges
 523
 524   Float_t xIntersect, yIntersect;
 525   Float_t x1, x2, y1, y2;
 526
 527   Float_t shiftX = GetShiftX();
 528   Float_t shiftY = GetShiftY();
 529
 530   Float_t xPoint = GetXPointOnCathode() + shiftX;
 531   Float_t yPoint = GetYPointOnCathode() + shiftY;
 532
 533   Float_t xemiss = GetXCoordOfEmission();
 534   Float_t yemiss = GetYCoordOfEmission();
 535
 536   Float_t phi = GetPhiPhotonInDRS();
 537   Float_t m = tan(phi);
 538
 539   Float_t x0 = xemiss + shiftX;
 540   Float_t y0 = yemiss + shiftY;
 541
 542   if(xPoint > x0)
 543     {
 544       x1 = x0;
 545       x2 = xPoint;
 546     }
 547   else
 548     {
 549       x2 = x0;
 550       x1 = xPoint;
 551     }
 552   if(yPoint > y0)
 553     {
 554       y1 = y0;
 555       y2 = yPoint;
 556     }
 557   else
 558     {
 559       y2 = y0;
 560       y1 = yPoint;
 561     }
 562   //
 563   xIntersect = fXmax;
 564   yIntersect = m*(xIntersect - x0) + y0;
 565   if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
 566     {
 567       SetIntersectionX(xIntersect);
 568       SetIntersectionY(yIntersect);
 569       return;
 570     }
 571   //
 572   xIntersect = fXmin;
 573   yIntersect = m*(xIntersect - x0) + y0;
 574   if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
 575     {
 576       SetIntersectionX(xIntersect);
 577       SetIntersectionY(yIntersect);
 578       return;
 579     }
 580   //
 581   yIntersect = fYmax;
 582   xIntersect = (yIntersect - y0)/m + x0;
 583   if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
 584     {
 585       SetIntersectionX(xIntersect);
 586       SetIntersectionY(yIntersect);
 587       return;
 588     }
 589   //
 590   yIntersect = fYmin;
 591   xIntersect = (yIntersect - y0)/m + x0;
 592   if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
 593     {
 594       SetIntersectionX(xIntersect);
 595       SetIntersectionY(yIntersect);
 596       return;
 597     }
 598
 599   cout << " sono fuori!!!!!!" << endl;
 600
 601 }
 602 //__________________________________________________________________________________________________
 603 Int_t AliRICHRecon::CheckDetectorAcceptance() const
 604 {
 605   // check for the acceptance
 606
 607   // crosses X -2.6 2.6 cm
 608   // crosses Y -1 1 cm
 609
 610   Float_t xcoord = GetDetectorWhereX();
 611   Float_t ycoord = GetDetectorWhereY();
 612
 613   if(xcoord > fXmax)
 614     {
 615       if(ycoord > fYmax) return 2;
 616       if(ycoord > fYmin && ycoord < fYmax) return 3;
 617       if(ycoord < fYmin) return 4;
 618     }
 619   if(xcoord < fXmin)
 620     {
 621       if(ycoord > fYmax) return 8;
 622       if(ycoord > fYmin && ycoord < fYmax) return 7;
 623       if(ycoord < fYmin) return 6;
 624     }
 625   if(xcoord > fXmin && xcoord < fXmax)
 626     {
 627       if(ycoord > fYmax) return 1;
 628       if(ycoord > fYmin && ycoord < fYmax) return 0;
 629       if(ycoord < fYmin) return 5;
 630     }
 631   return 999;
 632 }
 633 //__________________________________________________________________________________________________
 634 void AliRICHRecon::FindPhotonAnglesInDRS()
 635 {
 636   // Setup the rotation matrix of the track...
 637
 638   TRotation mtheta;
 639   TRotation mphi;
 640   TRotation minv;
 641   TRotation mrot;
 642
 643   Float_t trackTheta = GetTrackTheta();
 644   Float_t trackPhi = GetTrackPhi();
 645
 646   mtheta.RotateY(trackTheta);
 647   mphi.RotateZ(trackPhi);
 648
 649   mrot = mphi * mtheta;
 650   //  minv = mrot.Inverse();
 651
 652   TVector3 photonInRadiator(1,1,1);
 653
 654   Float_t thetaCerenkov = GetThetaPhotonInTRS();
 655   Float_t phiCerenkov   = GetPhiPhotonInTRS();
 656
 657   photonInRadiator.SetTheta(thetaCerenkov);
 658   photonInRadiator.SetPhi(phiCerenkov);
 659   photonInRadiator = mrot * photonInRadiator;
 660   Float_t theta = photonInRadiator.Theta();
 661   Float_t phi = photonInRadiator.Phi();
 662   SetThetaPhotonInDRS(theta);
 663   SetPhiPhotonInDRS(phi);
 664
 665 }
 666
 667 Float_t AliRICHRecon::FromEmissionToCathode()
 668 {
 669   // trace from emission point to cathode
 670
 671   Float_t nfreon, nquartz, ngas;
 672
 673   SetFreonRefractiveIndex();
 674   SetQuartzRefractiveIndex();
 675   SetGasRefractiveIndex();
 676
 677   nfreon  = GetFreonRefractiveIndex();
 678   nquartz = GetQuartzRefractiveIndex();
 679   ngas    = GetGasRefractiveIndex();
 680
 681   Float_t trackTheta = GetTrackTheta();
 682   Float_t trackPhi = GetTrackPhi();
 683   Float_t lengthOfEmissionPoint = GetEmissionPoint();
 684
 685   Float_t theta = GetThetaPhotonInDRS();
 686   Float_t phi   = GetPhiPhotonInDRS();
 687
 688 //   cout << " Theta " << Theta << " Phi " << Phi << endl;
 689
 690   Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi);
 691   Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi);
 692
 693   SetXCoordOfEmission(xemiss);
 694   SetYCoordOfEmission(yemiss);
 695
 696   Float_t thetaquar = SnellAngle(nfreon, nquartz, theta);
 697
 698   if(thetaquar == 999.)
 699     {
 700       SetXPointOnCathode(999.);
 701       SetYPointOnCathode(999.);
 702       return thetaquar;
 703     }
 704
 705   Float_t thetagap  = SnellAngle( nquartz, ngas, thetaquar);
 706
 707   if(thetagap == 999.)
 708     {
 709       SetXPointOnCathode(999.);
 710       SetYPointOnCathode(999.);
 711       return thetagap;
 712     }
 713
 714   Float_t xw = (fRadiatorWidth - lengthOfEmissionPoint)*cos(phi)*tan(theta);
 715   Float_t xq = fQuartzWidth*cos(phi)*tan(thetaquar);
 716   Float_t xg = fGapWidth*cos(phi)*tan(thetagap);
 717   Float_t yw = (fRadiatorWidth - lengthOfEmissionPoint)*sin(phi)*tan(theta);
 718   Float_t yq = fQuartzWidth*sin(phi)*tan(thetaquar);
 719   Float_t yg = fGapWidth*sin(phi)*tan(thetagap);
 720
 721
 722   Float_t xtot = xemiss + xw + xq + xg;
 723   Float_t ytot = yemiss + yw + yq + yg;
 724
 725   SetXPointOnCathode(xtot);
 726   SetYPointOnCathode(ytot);
 727
 728
 729   Float_t distanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2)
 730                                     +TMath::Power(fPhotonLimitY,2));
 731
 732   return distanceFromEntrance;
 733
 734 }
 735
 736
 737 void AliRICHRecon::FindPhiPoint()
 738 {
 739   //find phi of generated point
 740
 741   Float_t xtoentr = GetEntranceX();
 742   Float_t ytoentr = GetEntranceY();
 743
 744   Float_t trackTheta = GetTrackTheta();
 745   Float_t trackPhi = GetTrackPhi();
 746
 747   Float_t emissionPoint = GetEmissionPoint();
 748
 749   Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi);
 750   Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi);
 751   Float_t phi = atan2(argY,argX);
 752
 753   SetPhiPoint(phi);
 754
 755 }
 756
 757 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
 758 {
 759   // cerenkov angle from n and beta
 760
 761 // Compute the cerenkov angle
 762
 763   Float_t thetacer;
 764
 765   if((n*beta)<1.) {
 766     thetacer = 999.;
 767     //    cout << " warning in Cerenkoangle !!!!!! " << endl;
 768     return thetacer;
 769   }
 770
 771   thetacer = acos (1./(n*beta));
 772   return thetacer;
 773 }
 774
 775 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
 776 {
 777   // Snell law
 778
 779 // Compute the Snell angle
 780
 781   Float_t sinrefractangle;
 782   Float_t refractangle;
 783
 784   sinrefractangle = (n1/n2)*sin(theta1);
 785
 786   if(sinrefractangle>1.) {
 787     //    cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
 788     refractangle = 999.;
 789     return refractangle;
 790   }
 791
 792   refractangle = asin(sinrefractangle);
 793   return refractangle;
 794 }
 795
 796
 797 void AliRICHRecon::HoughResponse()
 798 {
 799   //Hough response
 800
 801 // Implement Hough response pat. rec. method
 802
 803   Float_t *hCSspace;
 804
 805   int           bin=0;
 806   int           bin1=0;
 807   int           bin2=0;
 808   int           i, j, k, nCorrBand;
 809   float         hcs[750],hcsw[750];
 810   float         angle, weight;
 811   float         lowerlimit,upperlimit;
 812
 813   float         etaPeak[100];
 814
 815   int           nBin;
 816
 817   float etaPeakPos  = -1;
 818
 819   Int_t   etaPeakCount = -1;
 820
 821   Float_t thetaCerenkov = 0.;
 822
 823   nBin = (int)(0.5+fThetaMax/(fDTheta));
 824   nCorrBand = (int)(0.5+ fWindowWidth/(2 * fDTheta));
 825
 826   memset ((void *)hcs, 0, fThetaBin*sizeof(float));
 827   memset ((void *)hcsw, 0, fThetaBin*sizeof(float));
 828
 829   Int_t nPhotons = GetPhotonsNumber();
 830
 831   Int_t weightFlag = 0;
 832
 833   for (k=0; k< nPhotons; k++) {
 834
 835     SetPhotonIndex(k);
 836
 837     angle = GetPhotonEta();
 838
 839     if(angle == -999.) continue;
 840
 841     if (angle>=fThetaMin && angle<= fThetaMax)
 842
 843       {
 844
 845         bin = (int)(0.5+angle/(fDTheta));
 846
 847         bin1= bin-nCorrBand;
 848         bin2= bin+nCorrBand;
 849
 850         // calculate weights
 851
 852         if(fIsWEIGHT)
 853           {
 854             lowerlimit = ((Float_t)bin1)*fDTheta + 0.5*fDTheta;
 855             SetThetaOfRing(lowerlimit);
 856             FindAreaAndPortionOfRing();
 857             Float_t area1 = GetAreaOfRing();
 858
 859             upperlimit = ((Float_t)bin2)*fDTheta + 0.5*fDTheta;
 860             SetThetaOfRing(upperlimit);
 861             FindAreaAndPortionOfRing();
 862             Float_t area2 = GetAreaOfRing();
 863
 864             //      cout << "lowerlimit" << lowerlimit << "upperlimit " << upperlimit << endl;
 865             Float_t diffarea = area2 - area1;
 866
 867             if(diffarea>0)
 868               {
 869                 weight = 1./(area2-area1);
 870               }
 871             else
 872               {
 873                 weightFlag = 1;
 874                 weight = 1.;
 875               }
 876
 877             //      cout <<" low "<< lowerlimit << " up " << upperlimit <<
 878             //        " area1 " << area1 << " area2 " << area2 << " weight " << weight << endl;
 879
 880           }
 881         else
 882           {
 883             weight = 1.;
 884           }
 885
 886         SetPhotonWeight(weight);
 887
 888         //      cout << "weight..." << weight << endl;
 889
 890
 891         if (bin1<0)    bin1=0;
 892         if (bin2>nBin) bin2=nBin;
 893
 894         for (j=bin1; j<bin2; j++)
 895           {
 896             hcs[j] += 1;
 897             hcsw[j] += weight;
 898           }
 899       }
 900   }
 901
 902
 903   if(weightFlag == 0)
 904     {
 905       hCSspace = hcsw;
 906     }
 907   else
 908     {
 909       hCSspace = hcs;
 910       //      cout << " probems with weight...normal procedure adopted " << endl;
 911     }
 912
 913   HoughFiltering(hCSspace);
 914
 915   for (bin=0; bin <nBin; bin++) {
 916     angle = (bin+0.5) * (fDTheta);
 917     if (hCSspace[bin] && hCSspace[bin] > etaPeakPos) {
 918       etaPeakCount = 0;
 919       etaPeakPos = hCSspace[bin];
 920       etaPeak[0]=angle;
 921     }
 922     else {
 923       if (hCSspace[bin] == etaPeakPos) {
 924         etaPeak[++etaPeakCount] = angle;
 925       }
 926     }
 927   }
 928
 929   for (i=0; i<etaPeakCount+1; i++) {
 930     thetaCerenkov += etaPeak[i];
 931   }
 932   if (etaPeakCount>=0) {
 933     thetaCerenkov /= etaPeakCount+1;
 934     fThetaPeakPos = etaPeakPos;
 935   }
 936
 937   SetThetaCerenkov(thetaCerenkov);
 938 }
 939
 940
 941 void AliRICHRecon::HoughFiltering(float hcs[])
 942 {
 943   // filter for Hough
 944
 945 // hough filtering
 946
 947    float hcsFilt[750];
 948    float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
 949    int nx, i, nxDx;
 950    int sizeHCS;
 951    int nBin;
 952
 953    nBin =  (int)(1+fThetaMax/fDTheta);
 954    sizeHCS = fThetaBin*sizeof(float);
 955
 956    memset ((void *)hcsFilt, 0, sizeHCS);
 957
 958    for (nx = 0; nx < nBin; nx++) {
 959       for (i = 0; i < 5; i++)   {
 960         nxDx = nx + (i-2);
 961         if (nxDx> -1 && nxDx<nBin)
 962              hcsFilt[nx] +=  hcs[nxDx] * k[i];
 963       }
 964    }
 965
 966    for (nx = 0; nx < nBin; nx++) {
 967      hcs[nx] = hcsFilt[nx];
 968    }
 969 }
 970
 971 void AliRICHRecon::FindWeightThetaCerenkov()
 972 {
 973   // manage with weight for photons
 974
 975   Float_t wei = 0.;
 976   Float_t weightThetaCerenkov = 0.;
 977
 978   Int_t nPhotons = GetPhotonsNumber();
 979   for(Int_t i=0;i<nPhotons;i++)
 980     {
 981       SetPhotonIndex(i);
 982
 983       if(GetPhotonFlag() == 2)
 984         {
 985           Float_t photonEta = GetPhotonEta();
 986           Float_t photonWeight = GetPhotonWeight();
 987           weightThetaCerenkov += photonEta*photonWeight;
 988           wei += photonWeight;
 989         }
 990     }
 991
 992   if(wei != 0.)
 993     {
 994       weightThetaCerenkov /= wei;
 995     }
 996   else
 997     {
 998       weightThetaCerenkov = 0.;
 999     }
1000
1001   SetThetaCerenkov(weightThetaCerenkov);
1002
1003   cout << " thetac weighted -> " << weightThetaCerenkov << endl;
1004 }
1005
1006
1007 void AliRICHRecon::FlagPhotons()
1008 {
1009   // flag photons
1010
1011   Int_t nPhotonHough = 0;
1012
1013   Float_t thetaCerenkov = GetThetaCerenkov();
1014   AliDebug(1,Form(" fThetaCerenkov %f ",thetaCerenkov));
1015
1016   Float_t thetaDist= thetaCerenkov - fThetaMin;
1017   Int_t steps = (Int_t)(thetaDist / fDTheta);
1018
1019   Float_t tmin = fThetaMin + (Float_t)(steps - 1)*fDTheta;
1020   Float_t tmax = fThetaMin + (Float_t)(steps)*fDTheta;
1021   Float_t tavg = 0.5*(tmin+tmax);
1022
1023   tmin = tavg - 0.5*fWindowWidth;
1024   tmax = tavg + 0.5*fWindowWidth;
1025
1026   //  Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
1027
1028   Int_t nPhotons = GetPhotonsNumber();
1029
1030   //  for(Int_t i=0;i<candidatePhotonsNumber;i++)
1031
1032   for(Int_t i=0;i<nPhotons;i++)
1033     {
1034       SetPhotonIndex(i);
1035
1036       Float_t photonEta = GetPhotonEta();
1037
1038       if(photonEta == -999.) continue;
1039
1040       if(photonEta >= tmin && photonEta <= tmax)
1041         {
1042           SetPhotonFlag(2);
1043           nPhotonHough++;
1044         }
1045     }
1046   SetHoughPhotons(nPhotonHough);
1047 }
1048