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