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