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 -10;//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(-1);
  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(-1);
 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 = 0;
 519   if (nPsiTotal>0)
 520     portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal));
 521
 522
 523   SetAreaOfRing(area);
 524   SetPortionOfRing(portionOfRing);
 525 }
 526 //__________________________________________________________________________________________________
 527 void AliRICHRecon::FindIntersectionWithDetector()
 528 {
 529   // find ring intersection with CsI edges
 530
 531   Float_t xIntersect, yIntersect;
 532   Float_t x1, x2, y1, y2;
 533
 534   Float_t shiftX = GetShiftX();
 535   Float_t shiftY = GetShiftY();
 536
 537   Float_t xPoint = GetXPointOnCathode() + shiftX;
 538   Float_t yPoint = GetYPointOnCathode() + shiftY;
 539
 540   Float_t xemiss = GetXCoordOfEmission();
 541   Float_t yemiss = GetYCoordOfEmission();
 542
 543   Float_t phi = GetPhiPhotonInDRS();
 544   Float_t m = tan(phi);
 545
 546   Float_t x0 = xemiss + shiftX;
 547   Float_t y0 = yemiss + shiftY;
 548
 549   if(xPoint > x0)
 550     {
 551       x1 = x0;
 552       x2 = xPoint;
 553     }
 554   else
 555     {
 556       x2 = x0;
 557       x1 = xPoint;
 558     }
 559   if(yPoint > y0)
 560     {
 561       y1 = y0;
 562       y2 = yPoint;
 563     }
 564   else
 565     {
 566       y2 = y0;
 567       y1 = yPoint;
 568     }
 569   //
 570   xIntersect = fXmax;
 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   xIntersect = fXmin;
 580   yIntersect = m*(xIntersect - x0) + y0;
 581   if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
 582     {
 583       SetIntersectionX(xIntersect);
 584       SetIntersectionY(yIntersect);
 585       return;
 586     }
 587   //
 588   yIntersect = fYmax;
 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   yIntersect = fYmin;
 598   xIntersect = (yIntersect - y0)/m + x0;
 599   if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
 600     {
 601       SetIntersectionX(xIntersect);
 602       SetIntersectionY(yIntersect);
 603       return;
 604     }
 605
 606   cout << " sono fuori!!!!!!" << endl;
 607
 608 }
 609 //__________________________________________________________________________________________________
 610 Int_t AliRICHRecon::CheckDetectorAcceptance() const
 611 {
 612   // check for the acceptance
 613
 614   // crosses X -2.6 2.6 cm
 615   // crosses Y -1 1 cm
 616
 617   Float_t xcoord = GetDetectorWhereX();
 618   Float_t ycoord = GetDetectorWhereY();
 619
 620   if(xcoord > fXmax)
 621     {
 622       if(ycoord > fYmax) return 2;
 623       if(ycoord > fYmin && ycoord < fYmax) return 3;
 624       if(ycoord < fYmin) return 4;
 625     }
 626   if(xcoord < fXmin)
 627     {
 628       if(ycoord > fYmax) return 8;
 629       if(ycoord > fYmin && ycoord < fYmax) return 7;
 630       if(ycoord < fYmin) return 6;
 631     }
 632   if(xcoord > fXmin && xcoord < fXmax)
 633     {
 634       if(ycoord > fYmax) return 1;
 635       if(ycoord > fYmin && ycoord < fYmax) return 0;
 636       if(ycoord < fYmin) return 5;
 637     }
 638   return 999;
 639 }
 640 //__________________________________________________________________________________________________
 641 void AliRICHRecon::FindPhotonAnglesInDRS()
 642 {
 643   // Setup the rotation matrix of the track...
 644
 645   TRotation mtheta;
 646   TRotation mphi;
 647   TRotation minv;
 648   TRotation mrot;
 649
 650   Float_t trackTheta = GetTrackTheta();
 651   Float_t trackPhi = GetTrackPhi();
 652
 653   mtheta.RotateY(trackTheta);
 654   mphi.RotateZ(trackPhi);
 655
 656   mrot = mphi * mtheta;
 657   //  minv = mrot.Inverse();
 658
 659   TVector3 photonInRadiator(1,1,1);
 660
 661   Float_t thetaCerenkov = GetThetaPhotonInTRS();
 662   Float_t phiCerenkov   = GetPhiPhotonInTRS();
 663
 664   photonInRadiator.SetTheta(thetaCerenkov);
 665   photonInRadiator.SetPhi(phiCerenkov);
 666   photonInRadiator = mrot * photonInRadiator;
 667   Float_t theta = photonInRadiator.Theta();
 668   Float_t phi = photonInRadiator.Phi();
 669   SetThetaPhotonInDRS(theta);
 670   SetPhiPhotonInDRS(phi);
 671
 672 }
 673 //__________________________________________________________________________________________________
 674 Float_t AliRICHRecon::FromEmissionToCathode()
 675 {
 676   // trace from emission point to cathode
 677
 678   Float_t nfreon, nquartz, ngas;
 679
 680   SetFreonRefractiveIndex();
 681   SetQuartzRefractiveIndex();
 682   SetGasRefractiveIndex();
 683
 684   nfreon  = GetFreonRefractiveIndex();
 685   nquartz = GetQuartzRefractiveIndex();
 686   ngas    = GetGasRefractiveIndex();
 687
 688   Float_t trackTheta = GetTrackTheta();
 689   Float_t trackPhi = GetTrackPhi();
 690   Float_t lengthOfEmissionPoint = GetEmissionPoint();
 691
 692   Float_t theta = GetThetaPhotonInDRS();
 693   Float_t phi   = GetPhiPhotonInDRS();
 694
 695 //   cout << " Theta " << Theta << " Phi " << Phi << endl;
 696
 697   Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi);
 698   Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi);
 699
 700   SetXCoordOfEmission(xemiss);
 701   SetYCoordOfEmission(yemiss);
 702
 703   Float_t thetaquar = SnellAngle(nfreon, nquartz, theta);
 704
 705   if(thetaquar == 999.)
 706     {
 707       SetXPointOnCathode(999.);
 708       SetYPointOnCathode(999.);
 709       return thetaquar;
 710     }
 711
 712   Float_t thetagap  = SnellAngle( nquartz, ngas, thetaquar);
 713
 714   if(thetagap == 999.)
 715     {
 716       SetXPointOnCathode(999.);
 717       SetYPointOnCathode(999.);
 718       return thetagap;
 719     }
 720
 721   Float_t xw = (fRadiatorWidth - lengthOfEmissionPoint)*cos(phi)*tan(theta);
 722   Float_t xq = fQuartzWidth*cos(phi)*tan(thetaquar);
 723   Float_t xg = fGapWidth*cos(phi)*tan(thetagap);
 724   Float_t yw = (fRadiatorWidth - lengthOfEmissionPoint)*sin(phi)*tan(theta);
 725   Float_t yq = fQuartzWidth*sin(phi)*tan(thetaquar);
 726   Float_t yg = fGapWidth*sin(phi)*tan(thetagap);
 727
 728
 729   Float_t xtot = xemiss + xw + xq + xg;
 730   Float_t ytot = yemiss + yw + yq + yg;
 731
 732   SetXPointOnCathode(xtot);
 733   SetYPointOnCathode(ytot);
 734
 735
 736   Float_t distanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2)
 737                                     +TMath::Power(fPhotonLimitY,2));
 738
 739   return distanceFromEntrance;
 740
 741 }
 742 //__________________________________________________________________________________________________
 743 void AliRICHRecon::FindPhiPoint()
 744 {
 745   //find phi of generated point
 746
 747   Float_t xtoentr = GetEntranceX();
 748   Float_t ytoentr = GetEntranceY();
 749
 750   Float_t trackTheta = GetTrackTheta();
 751   Float_t trackPhi = GetTrackPhi();
 752
 753   Float_t emissionPoint = GetEmissionPoint();
 754
 755   Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi);
 756   Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi);
 757   Float_t phi = atan2(argY,argX);
 758
 759   SetPhiPoint(phi);
 760
 761 }
 762 //__________________________________________________________________________________________________
 763 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
 764 {
 765   // cerenkov angle from n and beta
 766
 767 // Compute the cerenkov angle
 768
 769   Float_t thetacer;
 770
 771   if((n*beta)<1.) {
 772     thetacer = 999.;
 773     //    cout << " warning in Cerenkoangle !!!!!! " << endl;
 774     return thetacer;
 775   }
 776
 777   thetacer = acos (1./(n*beta));
 778   return thetacer;
 779 }
 780 //__________________________________________________________________________________________________
 781 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
 782 {
 783   // Snell law
 784
 785 // Compute the Snell angle
 786
 787   Float_t sinrefractangle;
 788   Float_t refractangle;
 789
 790   sinrefractangle = (n1/n2)*sin(theta1);
 791
 792   if(sinrefractangle>1.) {
 793     //    cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
 794     refractangle = 999.;
 795     return refractangle;
 796   }
 797
 798   refractangle = asin(sinrefractangle);
 799   return refractangle;
 800 }
 801 //__________________________________________________________________________________________________
 802 void AliRICHRecon::HoughResponse()
 803 {
 804   Int_t nChannels = (Int_t)(fThetaMax/fDTheta+0.5);
 805   TH1F *phots = new TH1F("phots","phots",nChannels,0.,fThetaMax);
 806   TH1F *photsw = new TH1F("photsw","photsw",nChannels,0.,fThetaMax);
 807   TH1F *resultw = new TH1F("resultw","resultw",nChannels,0.,fThetaMax);
 808   Int_t nPhotons = GetPhotonsNumber();
 809   Int_t nBin = (Int_t)(fThetaMax/fDTheta);
 810   Int_t nCorrBand = (Int_t)(fWindowWidth/(2*fDTheta));
 811   AliDebug(1,Form("Ring reconstruction for track with theta %f",GetTrackTheta()*TMath::RadToDeg()));
 812   for (Int_t kPhot=0; kPhot< nPhotons; kPhot++){
 813     SetPhotonIndex(kPhot);
 814     Double_t angle = GetPhotonEta();
 815     if(angle<0||angle>fThetaMax) continue;
 816     phots->Fill(angle);
 817     Int_t bin = (Int_t)(0.5+angle/(fDTheta));
 818     Double_t weight=1.;
 819     if(fIsWEIGHT){
 820       Double_t lowerlimit = ((Float_t)bin)*fDTheta - 0.5*fDTheta;
 821       SetThetaOfRing(lowerlimit);
 822       FindAreaAndPortionOfRing();
 823       Float_t area1 = GetAreaOfRing();
 824       Double_t upperlimit = ((Float_t)bin)*fDTheta + 0.5*fDTheta;
 825       SetThetaOfRing(upperlimit);
 826       FindAreaAndPortionOfRing();
 827       Float_t area2 = GetAreaOfRing();
 828       AliDebug(1,Form("lowerlimit %f  area %f ; upperlimit %f area %f",lowerlimit,area1,upperlimit,area2));
 829       Float_t diffarea = area2 - area1;
 830       if(diffarea>0){weight = 1./(area2-area1);}else{weight = 1.;}
 831     }
 832     AliDebug(1,Form("Calculated weight %f",weight));
 833     photsw->Fill(angle,weight);
 834     SetPhotonWeight(weight);
 835   }
 836   for (Int_t i=1; i<=nBin;i++){
 837     Int_t bin1= i-nCorrBand;
 838     Int_t bin2= i+nCorrBand;
 839     if(bin1<1) bin1=1;
 840     if(bin2>nBin)bin2=nBin;
 841     Double_t sumPhots=phots->Integral(bin1,bin2);
 842     if(sumPhots<fMinNumPhots) continue; // cut on minimum n. of photons per ring
 843     Double_t sumPhotsw=photsw->Integral(bin1,bin2);
 844     resultw->Fill((Float_t)((i+0.5)*fDTheta),sumPhotsw);
 845 }
 846 // evaluate the "BEST" thetacerenkov....
 847   Float_t *pVec = resultw->GetArray();
 848   Int_t locMax = TMath::LocMax(nBin,pVec);
 849   SetThetaCerenkov((Double_t)(locMax*fDTheta+0.5*fDTheta));
 850
 851 // Reset and delete objects
 852   phots->Delete();photsw->Delete();resultw->Delete();
 853 }//HoughResponse
 854 //__________________________________________________________________________________________________
 855 void AliRICHRecon::FindThetaCerenkov()
 856 {
 857   // manage with weight for photons
 858
 859   Float_t wei = 0.;
 860   Float_t weightThetaCerenkov = 0.;
 861
 862   Int_t nPhotons = GetPhotonsNumber();
 863   Double_t etaMin=9999.,etaMax=0.;
 864   for(Int_t i=0;i<nPhotons;i++)
 865     {
 866       SetPhotonIndex(i);
 867
 868       if(GetPhotonFlag() == 2)
 869         {
 870           Float_t photonEta = GetPhotonEta();
 871           if(photonEta<etaMin) etaMin=photonEta;
 872           if(photonEta>etaMax) etaMax=photonEta;
 873           Float_t photonWeight = GetPhotonWeight();
 874           weightThetaCerenkov += photonEta*photonWeight;
 875           wei += photonWeight;
 876         }
 877     }
 878
 879   if(wei != 0.) weightThetaCerenkov /= wei; else weightThetaCerenkov = 0.;
 880   SetThetaCerenkov(weightThetaCerenkov);
 881
 882   // estimate of the n. of bkg photons
 883   SetThetaOfRing(etaMin); FindAreaAndPortionOfRing(); Double_t internalArea = GetAreaOfRing();
 884   SetThetaOfRing(etaMax); FindAreaAndPortionOfRing(); Double_t externalArea = GetAreaOfRing();
 885
 886   Double_t effArea = (AliRICHParam::PcSizeX()-AliRICHParam::DeadZone())*(AliRICHParam::PcSizeY()-2*AliRICHParam::DeadZone());
 887   Double_t nPhotBKG = (externalArea-internalArea)/effArea*fpClusters->GetEntries();
 888   if(nPhotBKG<0) nPhotBKG=0; //just protection from funny angles...
 889   SetPhotBKG(nPhotBKG);
 890   //
 891
 892   AliDebug(1,Form(" thetac weighted -> %f",weightThetaCerenkov));
 893 }
 894 //__________________________________________________________________________________________________
 895 void AliRICHRecon::FlagPhotons()
 896 {
 897   // flag photons
 898
 899   Int_t nPhotonHough = 0;
 900
 901   Float_t thetaCerenkov = GetThetaCerenkov();
 902   AliDebug(1,Form(" fThetaCerenkov %f ",thetaCerenkov));
 903
 904   Float_t thetaDist= thetaCerenkov - fThetaMin;
 905   Int_t steps = (Int_t)(thetaDist / fDTheta);
 906
 907   Float_t tmin = fThetaMin + (Float_t)(steps - 1)*fDTheta;
 908   Float_t tmax = fThetaMin + (Float_t)(steps)*fDTheta;
 909   Float_t tavg = 0.5*(tmin+tmax);
 910
 911   tmin = tavg - 0.5*fWindowWidth;
 912   tmax = tavg + 0.5*fWindowWidth;
 913
 914   //  Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
 915
 916   Int_t nPhotons = GetPhotonsNumber();
 917
 918   //  for(Int_t i=0;i<candidatePhotonsNumber;i++)
 919
 920 //  TH1F* h1_phots = (TH1F*)gROOT->FindObjectAny("h1_phots");
 921
 922 //  cout << "h1_phots " << h1_phots << endl;
 923
 924   for(Int_t i=0;i<nPhotons;i++)
 925     {
 926       SetPhotonIndex(i);
 927
 928       Float_t photonEta = GetPhotonEta();
 929
 930       if(photonEta == -999.) continue;
 931
 932       if(photonEta >= tmin && photonEta <= tmax)
 933         {
 934           SetPhotonFlag(2);
 935           nPhotonHough++;
 936 //          if(h1_phots)h1_phots->Fill(photonEta);
 937         }
 938     }
 939   SetHoughPhotons(nPhotonHough);
 940 }