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