1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //////////////////////////////////////////////////////////////////////////
20 // RICH class to perfom pattern recognition based on Hough transfrom //
21 // for single chamber //
22 //////////////////////////////////////////////////////////////////////////
24 #include "AliRICHRecon.h" //class header
26 #include <TRotation.h>
30 #include "AliRICHCluster.h" //ThetaCerenkov()
31 #include "AliRICHParam.h"
32 #include "AliRICHHelix.h" //ThetaCerenkov()
35 #define NPointsOfRing 201
37 //__________________________________________________________________________________________________
38 AliRICHRecon::AliRICHRecon():
39 TTask ("RichRec","RichPat"),
41 fPhotonIndex(0), // should we use -1?
42 fFittedHoughPhotons(0),
64 fLengthEmissionPoint(0),
67 fDistanceFromCluster(999),
68 fCerenkovAnglePad(999),
69 fThetaPhotonCerenkov(0),
76 fMassHypotesis(0.139567),
84 fFittedThetaCerenkov(0),
88 fIsBACKGROUND(kFALSE),
94 fWeightThetaCerenkov(0),
99 for (Int_t i=0; i<3000; i++) {
102 fPhotonEta[i] = -999;
103 fPhotonWeight[i] = 0.0;
106 fWeightPhotons[i] = 0;
109 //__________________________________________________________________________________________________
110 Double_t AliRICHRecon::ThetaCerenkov(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t &nphot)
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.
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;
124 SetThetaCerenkov(-1);
127 // Photon Flag: Flag = 0 initial set; Flag = 1 good candidate (charge compatible with photon); Flag = 2 photon used for the ring;
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);
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
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);
148 SetPhotonsNumber(pClusters->GetEntries());
150 if((nphot=FlagPhotons(HoughResponse()))<1) return -11; //flag photons according to individual theta ckov with respect to most probable track theta ckov
153 // Float_t thetaCerenkov = GetThetaCerenkov();
154 // SetThetaOfRing(thetaCerenkov);
155 // FindAreaAndPortionOfRing();
157 // Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
158 // SetHoughPhotonsNorm(nPhotonHoughNorm);
160 // Calculate the area where the photon are accepted...
162 Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth;
163 SetThetaOfRing(thetaInternal);
164 FindAreaAndPortionOfRing();
165 Float_t internalArea = GetAreaOfRing();
167 Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth;
168 SetThetaOfRing(thetaExternal);
169 FindAreaAndPortionOfRing();
170 Float_t externalArea = GetAreaOfRing();
172 Float_t houghArea = externalArea - internalArea;
174 SetHoughArea(houghArea);
176 FindThetaCerenkov(pClusters->GetEntries());
177 return GetThetaCerenkov();
179 //__________________________________________________________________________________________________
180 Int_t AliRICHRecon::PhotonInBand()
182 // Define valid band for photon candidates. For that photons with ThetaMin and ThetaMax are traced up to photcathode
188 Float_t xtoentr = GetEntranceX();
189 Float_t ytoentr = GetEntranceY();
194 Float_t phpad = GetPhiPoint();
198 SetEmissionPoint(AliRICHParam::RadThick() -0.0001);
200 nfreon = AliRICHParam::Instance()->IdxC6F14(AliRICHParam::Instance()->EckovMin());
203 AliDebug(1,Form("thetacer in photoninband min %f",thetacer));
205 FindThetaAtQuartz(thetacer);
207 if(thetacer == 999. || GetThetaAtQuartz() == 999.) {
209 SetXInnerRing(-999.);
210 SetYInnerRing(-999.);
211 SetRadiusInnerRing(-999.);
215 SetThetaPhotonInDRS(GetThetaAtQuartz());
216 SetPhiPhotonInDRS(phpad);
218 innerRadius = FromEmissionToCathode();
219 if(innerRadius == 999.) innerRadius = -999.;
221 SetXInnerRing(GetXPointOnCathode());
222 SetYInnerRing(GetYPointOnCathode());
223 SetRadiusInnerRing(innerRadius);
227 SetEmissionPoint(0.);
228 // SetMassHypotesis(0.139567);
230 nfreon = AliRICHParam::Instance()->IdxC6F14(AliRICHParam::Instance()->EckovMax());
232 thetacer = Cerenkovangle(nfreon,1);
236 AliDebug(1,Form("thetacer in photoninband max %f",thetacer));
238 FindThetaAtQuartz(thetacer);
240 if(thetacer == 999. || GetThetaAtQuartz() == 999.)
245 SetRadiusOuterRing(999.);
249 SetThetaPhotonInDRS(GetThetaAtQuartz());
250 SetPhiPhotonInDRS(phpad);
252 outerRadius = FromEmissionToCathode();
253 // cout << " outerRadius " << outerRadius << endl;
254 SetXOuterRing(GetXPointOnCathode());
255 SetYOuterRing(GetYPointOnCathode());
256 SetRadiusOuterRing(outerRadius);
259 Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2));
261 AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius));
263 if(padradius>=innerRadius && padradius<=outerRadius) return 1;
266 //__________________________________________________________________________________________________
267 void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
269 // find the theta at the quartz plate
271 if(thetaCerenkov == 999.) { SetThetaAtQuartz(999.); return; }
273 Float_t thetaAtQuartz = 999.;
275 Float_t trackTheta = GetTrackTheta();
277 if(trackTheta == 0) {
278 thetaAtQuartz = thetaCerenkov;
279 SetThetaAtQuartz(thetaAtQuartz);
283 Float_t trackPhi = GetTrackPhi();
284 Float_t phiPoint = GetPhiPoint();
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;
291 Double_t underSqrt = 1 + b*b - c*c;
294 SetThetaAtQuartz(999.);
298 Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c);
299 Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c);
301 Double_t thetaSol1 = 2*TMath::ATan(sol1);
302 Double_t thetaSol2 = 2*TMath::ATan(sol2);
304 if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1;
305 if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2;
307 // AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz));
309 SetThetaAtQuartz(thetaAtQuartz);
311 //__________________________________________________________________________________________________
312 void AliRICHRecon::FindThetaPhotonCerenkov()
314 //find theta cerenkov of ring
316 Float_t thetaCerMin = 0.;
317 Float_t thetaCerMax = 0.75;
318 Float_t thetaCerMean;
320 Float_t radiusMin, radiusMax, radiusMean;
321 Int_t nIteration = 0;
323 const Float_t kTollerance = 0.05;
326 Float_t phiPoint = GetPhiPoint();
328 SetEmissionPoint(AliRICHParam::RadThick()/2);
330 Float_t xPoint = GetEntranceX();
331 Float_t yPoint = GetEntranceY();
332 Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint);
334 // AliDebug(1,Form(" DistPoint %f ",distPoint));
336 // Star minimization...
340 FindThetaAtQuartz(thetaCerMin);
342 if(GetThetaAtQuartz() == 999.)
348 SetThetaPhotonInDRS(GetThetaAtQuartz());
349 SetPhiPhotonInDRS(phiPoint);
351 radiusMin = FromEmissionToCathode();
356 FindThetaAtQuartz(thetaCerMax);
357 if(GetThetaAtQuartz() == 999.)
363 SetThetaPhotonInDRS(GetThetaAtQuartz());
364 SetPhiPhotonInDRS(phiPoint);
366 radiusMax = FromEmissionToCathode();
370 thetaCerMean = (thetaCerMax + thetaCerMin)/2;
372 FindThetaAtQuartz(thetaCerMean);
373 if(GetThetaAtQuartz() == 999.)
379 SetThetaPhotonInDRS(GetThetaAtQuartz());
380 SetPhiPhotonInDRS(phiPoint);
382 radiusMean = FromEmissionToCathode();
385 // AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax));
387 while (TMath::Abs(radiusMean-distPoint) > kTollerance)
390 if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean;
391 if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) {
393 thetaCerMin = thetaCerMean;
395 FindThetaAtQuartz(thetaCerMin);
396 SetThetaPhotonInDRS(GetThetaAtQuartz());
397 SetPhiPhotonInDRS(phiPoint);
399 radiusMin =FromEmissionToCathode();
402 thetaCerMean = (thetaCerMax + thetaCerMin)/2;
404 FindThetaAtQuartz(thetaCerMean);
405 SetThetaPhotonInDRS(GetThetaAtQuartz());
406 SetPhiPhotonInDRS(phiPoint);
408 radiusMean = FromEmissionToCathode();
412 // AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration));
413 SetThetaPhotonCerenkov(999.);
418 // AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean));
419 SetThetaPhotonCerenkov(thetaCerMean);
422 //__________________________________________________________________________________________________
423 void AliRICHRecon::FindAreaAndPortionOfRing()
425 //find fraction of the ring accepted by the RICH
427 Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing];
429 // Float_t xtoentr = GetEntranceX();
430 // Float_t ytoentr = GetEntranceY();
431 Float_t shiftX = GetShiftX();
432 Float_t shiftY = GetShiftY();
434 Float_t xemiss = GetXCoordOfEmission();
435 Float_t yemiss = GetYCoordOfEmission();
437 Float_t x0 = xemiss + shiftX;
438 Float_t y0 = yemiss + shiftY;
442 SetEmissionPoint(AliRICHParam::RadThick()/2.);
444 Float_t theta = GetThetaOfRing();
447 Int_t nPsiAccepted = 0;
450 for(Int_t i=0;i<NPointsOfRing-1;i++){
451 Float_t psi = 2*TMath::Pi()*i/NPointsOfRing;
453 SetThetaPhotonInTRS(theta);
454 SetPhiPhotonInTRS(psi);
455 FindPhotonAnglesInDRS();
457 Float_t radius = FromEmissionToCathode();
458 if (radius == 999.) continue;
462 Float_t xPointRing = GetXPointOnCathode() + shiftX;
463 Float_t yPointRing = GetYPointOnCathode() + shiftY;
465 SetDetectorWhereX(xPointRing);
466 SetDetectorWhereY(yPointRing);
468 Int_t zone = CheckDetectorAcceptance();
470 AliDebug(1,Form("acceptance to detector zone -> %d",zone));
473 FindIntersectionWithDetector();
474 xPoint[nPoints] = GetIntersectionX(); yPoint[nPoints] = GetIntersectionY();
476 xPoint[nPoints] = xPointRing; yPoint[nPoints] = yPointRing;
482 xPoint[nPoints] = xPoint[0]; yPoint[nPoints] = yPoint[0];
488 for (Int_t i = 0; i < nPoints; i++)
490 area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0));
495 Float_t portionOfRing = 0;
497 portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal));
501 SetPortionOfRing(portionOfRing);
502 }//FindAreaAndPortionOfRing()
503 //__________________________________________________________________________________________________
504 void AliRICHRecon::FindIntersectionWithDetector()
506 // find ring intersection with CsI edges
508 Float_t xIntersect, yIntersect;
509 Float_t x1, x2, y1, y2;
511 Float_t shiftX = GetShiftX();
512 Float_t shiftY = GetShiftY();
514 Float_t xPoint = GetXPointOnCathode() + shiftX;
515 Float_t yPoint = GetYPointOnCathode() + shiftY;
517 Float_t xemiss = GetXCoordOfEmission();
518 Float_t yemiss = GetYCoordOfEmission();
520 Float_t phi = GetPhiPhotonInDRS();
521 Float_t m = tan(phi);
523 Float_t x0 = xemiss + shiftX;
524 Float_t y0 = yemiss + shiftY;
547 xIntersect = AliRICHParam::PcSizeX();
548 yIntersect = m*(xIntersect - x0) + y0;
549 if (yIntersect >= 0 && yIntersect <= AliRICHParam::PcSizeY() && xIntersect >= x1 && xIntersect <= x2)
551 SetIntersectionX(xIntersect);
552 SetIntersectionY(yIntersect);
557 yIntersect = m*(xIntersect - x0) + y0;
558 if (yIntersect >= 0 && yIntersect <= AliRICHParam::PcSizeY() && xIntersect >= x1 && xIntersect <= x2)
560 SetIntersectionX(xIntersect);
561 SetIntersectionY(yIntersect);
565 yIntersect = AliRICHParam::PcSizeY();
566 xIntersect = (yIntersect - y0)/m + x0;
567 if (xIntersect >= 0 && xIntersect <= AliRICHParam::PcSizeX() && yIntersect >= y1 && yIntersect <= y2)
569 SetIntersectionX(xIntersect);
570 SetIntersectionY(yIntersect);
575 xIntersect = (yIntersect - y0)/m + x0;
576 if (xIntersect >= 0 && xIntersect <= AliRICHParam::PcSizeX() && yIntersect >= y1 && yIntersect <= y2)
578 SetIntersectionX(xIntersect);
579 SetIntersectionY(yIntersect);
584 //__________________________________________________________________________________________________
585 Int_t AliRICHRecon::CheckDetectorAcceptance() const
587 // check for the acceptance
589 // crosses X -2.6 2.6 cm
592 Float_t xcoord = GetDetectorWhereX();
593 Float_t ycoord = GetDetectorWhereY();
595 if(xcoord > AliRICHParam::PcSizeX())
597 if(ycoord > AliRICHParam::PcSizeY()) return 2;
598 if(ycoord > 0 && ycoord < AliRICHParam::PcSizeY()) return 3;
599 if(ycoord < 0) return 4;
603 if(ycoord > AliRICHParam::PcSizeY()) return 8;
604 if(ycoord > 0 && ycoord < AliRICHParam::PcSizeY()) return 7;
605 if(ycoord < 0) return 6;
607 if(xcoord > 0 && xcoord < AliRICHParam::PcSizeX())
609 if(ycoord > AliRICHParam::PcSizeY()) return 1;
610 if(ycoord > 0 && ycoord < AliRICHParam::PcSizeY()) return 0;
611 if(ycoord < 0) return 5;
615 //__________________________________________________________________________________________________
616 void AliRICHRecon::FindPhotonAnglesInDRS()
618 // Setup the rotation matrix of the track...
625 Float_t trackTheta = GetTrackTheta();
626 Float_t trackPhi = GetTrackPhi();
628 mtheta.RotateY(trackTheta);
629 mphi.RotateZ(trackPhi);
631 mrot = mphi * mtheta;
632 // minv = mrot.Inverse();
634 TVector3 photonInRadiator(1,1,1);
636 Float_t thetaCerenkov = GetThetaPhotonInTRS();
637 Float_t phiCerenkov = GetPhiPhotonInTRS();
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);
648 //__________________________________________________________________________________________________
649 Float_t AliRICHRecon::FromEmissionToCathode()
651 // Trace current photon from emission point somewhere in radiator to photocathode
655 Float_t nfreon, nquartz, ngas;
657 nfreon = AliRICHParam::Instance()->IdxC6F14(AliRICHParam::Instance()->EckovMean());
658 nquartz = AliRICHParam::Instance()->IdxSiO2(AliRICHParam::Instance()->EckovMean());
659 ngas = AliRICHParam::Instance()->IdxCH4(AliRICHParam::Instance()->EckovMean());
661 Float_t trackTheta = GetTrackTheta();
662 Float_t trackPhi = GetTrackPhi();
663 Float_t lengthOfEmissionPoint = GetEmissionPoint();
665 Float_t theta = GetThetaPhotonInDRS();
666 Float_t phi = GetPhiPhotonInDRS();
668 Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi);
669 Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi);
671 SetXCoordOfEmission(xemiss);
672 SetYCoordOfEmission(yemiss);
674 Float_t thetaquar = SnellAngle(nfreon, nquartz, theta);
676 if(thetaquar == 999.)
678 SetXPointOnCathode(999.);
679 SetYPointOnCathode(999.);
683 Float_t thetagap = SnellAngle( nquartz, ngas, thetaquar);
687 SetXPointOnCathode(999.);
688 SetYPointOnCathode(999.);
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);
700 Float_t xtot = xemiss + xw + xq + xg;
701 Float_t ytot = yemiss + yw + yq + yg;
703 SetXPointOnCathode(xtot);
704 SetYPointOnCathode(ytot);
707 Float_t distanceFromEntrance = TMath::Sqrt(TMath::Power(fPhotonLimitX,2)+TMath::Power(fPhotonLimitY,2));
709 return distanceFromEntrance;
712 //__________________________________________________________________________________________________
713 void AliRICHRecon::FindPhiPoint()
715 //find phi of generated point
717 Float_t xtoentr = GetEntranceX();
718 Float_t ytoentr = GetEntranceY();
720 Float_t trackTheta = GetTrackTheta();
721 Float_t trackPhi = GetTrackPhi();
723 Float_t emissionPoint = GetEmissionPoint();
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);
732 //__________________________________________________________________________________________________
733 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
735 // cerenkov angle from n and beta
737 // Compute the cerenkov angle
743 // cout << " warning in Cerenkoangle !!!!!! " << endl;
747 thetacer = acos (1./(n*beta));
750 //__________________________________________________________________________________________________
751 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
755 // Compute the Snell angle
757 Float_t sinrefractangle;
758 Float_t refractangle;
760 sinrefractangle = (n1/n2)*sin(theta1);
762 if(sinrefractangle>1.) {
763 // cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
768 refractangle = asin(sinrefractangle);
771 //__________________________________________________________________________________________________
772 Double_t AliRICHRecon::HoughResponse()
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;
789 Int_t bin = (Int_t)(0.5+angle/(fDTheta));
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.;}
804 AliDebug(1,Form("Calculated weight %f",weight));
805 photsw->Fill(angle,weight);
806 SetPhotonWeight(weight);
808 for (Int_t i=1; i<=nBin;i++){
809 Int_t bin1= i-nCorrBand;
810 Int_t bin2= i+nCorrBand;
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);
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
823 return (Double_t)(locMax*fDTheta+0.5*fDTheta); //final most probable track theta ckov
825 //__________________________________________________________________________________________________
826 void AliRICHRecon::FindThetaCerenkov(Int_t iNclus)
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
834 Float_t weightThetaCerenkov = 0.;
836 Double_t etaMin=9999.,etaMax=0.;
837 Double_t sigma2 = 0; //to collect error squared for this ring
839 for(Int_t i=0;i<GetPhotonsNumber();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;
848 //here comes sigma of the reconstructed ring
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();
855 Double_t beta = 1./(TMath::Cos(GetPhotonEta())*AliRICHParam::Instance()->IdxC6F14(AliRICHParam::EckovMean()));
856 sigma2 += 1./AliRICHParam::SigmaSinglePhotonFormula(GetPhotonEta(),GetPhiPoint(),GetTrackTheta(),phiref,beta);
860 if(sigma2>0) SetRingSigma2(1./sigma2);
861 else SetRingSigma2(1e10);
863 if(wei != 0.) weightThetaCerenkov /= wei; else weightThetaCerenkov = 0.;
864 SetThetaCerenkov(weightThetaCerenkov);
866 // estimate of the n. of bkg photons
867 SetThetaOfRing(etaMin); FindAreaAndPortionOfRing(); Double_t internalArea = GetAreaOfRing();
868 SetThetaOfRing(etaMax); FindAreaAndPortionOfRing(); Double_t externalArea = GetAreaOfRing();
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);
875 AliDebug(1,Form(" thetac weighted -> %f",weightThetaCerenkov));
876 }//FindThetaCerenkov()
877 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
878 Int_t AliRICHRecon::FlagPhotons(Double_t thetaCkovHough)
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
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
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);
890 tmin = tavg - 0.5*fWindowWidth; tmax = tavg + 0.5*fWindowWidth;
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) {