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