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