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