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 //
22 //////////////////////////////////////////////////////////////////////////
24 #include <Riostream.h>
26 #include <TRotation.h>
31 #include "AliRICHParam.h"
32 #include "AliRICHRecon.h"
33 #include "AliRICHHelix.h"
36 #define NPointsOfRing 201
38 //__________________________________________________________________________________________________
39 AliRICHRecon::AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId)
40 :TTask("RichRec","RichPat")
43 SetFreonScaleFactor(1);
44 // fIsWEIGHT = kFALSE;
46 fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75;
47 fDTheta = 0.001; fWindowWidth = 0.045;
49 fRadiatorWidth = AliRICHParam::Zfreon();
50 fQuartzWidth = AliRICHParam::Zwin();
51 fGapWidth = AliRICHParam::Freon2Pc() - fRadiatorWidth - fQuartzWidth;
53 fXmax = AliRICHParam::PcSizeX();
55 fYmax = AliRICHParam::PcSizeY();
56 SetTrackTheta(pHelix->Ploc().Theta());
57 SetTrackPhi(pHelix->Ploc().Phi());
59 SetShiftX(pHelix->PosRad().X());
60 SetShiftY(pHelix->PosRad().Y());
61 fpClusters = pClusters;
63 //__________________________________________________________________________________________________
64 Double_t AliRICHRecon::ThetaCerenkov()
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
69 if(fpClusters->GetEntries()==0) return -1;//no clusters at all for a given track
70 Bool_t kPatRec = kFALSE;
72 AliDebug(1,Form("---Track Parameters--- Theta: %f , Phi: %f ",GetTrackTheta()*TMath::RadToDeg(),GetTrackPhi()*TMath::RadToDeg()));
74 Int_t candidatePhotons = 0;
76 SetThetaCerenkov(999.);
78 SetHoughPhotonsNorm(0);
80 for (Int_t j=0; j < fpClusters->GetEntries(); j++){//clusters loop
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);
91 // Int_t photonStatus = PhotonInBand();
92 // if(photonStatus == 0) continue;
94 FindThetaPhotonCerenkov();
95 Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
96 AliDebug(1,Form("THETA CERENKOV ---> %f",thetaPhotonCerenkov));
97 SetPhotonEta(thetaPhotonCerenkov);
101 if(candidatePhotons >= 1) kPatRec = kTRUE;
103 if(!kPatRec) return -1;
105 SetPhotonsNumber(fpClusters->GetEntries());
112 Int_t nPhotonHough = GetHoughPhotons();
116 SetThetaCerenkov(999.);
117 SetHoughPhotonsNorm(0.);
123 AliDebug(1,Form("Number of clusters accepted ---> %i",nPhotonHough));
125 // Float_t thetaCerenkov = GetThetaCerenkov();
126 // SetThetaOfRing(thetaCerenkov);
127 // FindAreaAndPortionOfRing();
129 // Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
130 // SetHoughPhotonsNorm(nPhotonHoughNorm);
132 // Calculate the area where the photon are accepted...
134 Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth;
135 SetThetaOfRing(thetaInternal);
136 FindAreaAndPortionOfRing();
137 Float_t internalArea = GetAreaOfRing();
139 Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth;
140 SetThetaOfRing(thetaExternal);
141 FindAreaAndPortionOfRing();
142 Float_t externalArea = GetAreaOfRing();
144 Float_t houghArea = externalArea - internalArea;
146 SetHoughArea(houghArea);
148 return GetThetaCerenkov();
151 //__________________________________________________________________________________________________
152 void AliRICHRecon::FindEmissionPoint()
154 //estimate the emission point in radiator
156 // Find emission point
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;
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;
168 SetEmissionPoint(emissionPoint);
169 SetEmissionPoint(fRadiatorWidth/2); // tune the emission point
171 //__________________________________________________________________________________________________
172 Int_t AliRICHRecon::PhotonInBand()
174 //search band fro photon candidates
176 // Float_t massOfParticle;
181 Float_t xtoentr = GetEntranceX();
182 Float_t ytoentr = GetEntranceY();
187 Float_t phpad = GetPhiPoint();
191 SetPhotonEnergy(5.6);
192 SetEmissionPoint(fRadiatorWidth -0.0001);
193 SetFreonRefractiveIndex();
195 nfreon = GetFreonRefractiveIndex();
198 AliDebug(1,Form("thetacer in photoninband min %f",thetacer));
200 FindThetaAtQuartz(thetacer);
202 if(thetacer == 999. || GetThetaAtQuartz() == 999.)
205 SetXInnerRing(-999.);
206 SetYInnerRing(-999.);
207 SetRadiusInnerRing(-999.);
211 SetThetaPhotonInDRS(GetThetaAtQuartz());
212 SetPhiPhotonInDRS(phpad);
214 innerRadius = FromEmissionToCathode();
215 if(innerRadius == 999.) innerRadius = -999.;
217 SetXInnerRing(GetXPointOnCathode());
218 SetYInnerRing(GetYPointOnCathode());
219 SetRadiusInnerRing(innerRadius);
223 SetPhotonEnergy(7.7);
224 SetEmissionPoint(0.);
225 // SetMassHypotesis(0.139567);
226 SetFreonRefractiveIndex();
228 nfreon = GetFreonRefractiveIndex();
230 thetacer = Cerenkovangle(nfreon,1);
234 AliDebug(1,Form("thetacer in photoninband max %f",thetacer));
236 FindThetaAtQuartz(thetacer);
238 if(thetacer == 999. || GetThetaAtQuartz() == 999.)
243 SetRadiusOuterRing(999.);
247 SetThetaPhotonInDRS(GetThetaAtQuartz());
248 SetPhiPhotonInDRS(phpad);
250 outerRadius = FromEmissionToCathode();
251 // cout << " outerRadius " << outerRadius << endl;
252 SetXOuterRing(GetXPointOnCathode());
253 SetYOuterRing(GetYPointOnCathode());
254 SetRadiusOuterRing(outerRadius);
257 Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2));
259 AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius));
261 if(padradius>=innerRadius && padradius<=outerRadius) return 1;
264 //__________________________________________________________________________________________________
265 void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
267 //find the theta at the quartz plate
269 if(thetaCerenkov == 999.)
271 SetThetaAtQuartz(999.);
275 Float_t thetaAtQuartz = 999.;
277 Float_t trackTheta = GetTrackTheta();
279 if(trackTheta == 0) {
281 thetaAtQuartz = thetaCerenkov;
282 SetThetaAtQuartz(thetaAtQuartz);
286 Float_t trackPhi = GetTrackPhi();
287 Float_t phiPoint = GetPhiPoint();
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;
298 Double_t underSqrt = 1 + b*b - c*c;
301 SetThetaAtQuartz(999.);
305 Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c);
306 Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c);
308 Double_t thetaSol1 = 2*TMath::ATan(sol1);
309 Double_t thetaSol2 = 2*TMath::ATan(sol2);
311 if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1;
312 if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2;
314 // AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz));
316 SetThetaAtQuartz(thetaAtQuartz);
318 //__________________________________________________________________________________________________
319 void AliRICHRecon::FindThetaPhotonCerenkov()
321 //find theta cerenkov of ring
323 Float_t thetaCerMin = 0.;
324 Float_t thetaCerMax = 0.75;
325 Float_t thetaCerMean;
327 Float_t radiusMin, radiusMax, radiusMean;
328 Int_t nIteration = 0;
330 const Float_t kTollerance = 0.05;
333 Float_t phiPoint = GetPhiPoint();
335 SetPhotonEnergy(AliRICHParam::MeanCkovEnergy());
336 SetEmissionPoint(fRadiatorWidth/2);
338 Float_t xPoint = GetEntranceX();
339 Float_t yPoint = GetEntranceY();
340 Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint);
342 // AliDebug(1,Form(" DistPoint %f ",distPoint));
344 // Star minimization...
348 FindThetaAtQuartz(thetaCerMin);
350 if(GetThetaAtQuartz() == 999.)
356 SetThetaPhotonInDRS(GetThetaAtQuartz());
357 SetPhiPhotonInDRS(phiPoint);
359 radiusMin = FromEmissionToCathode();
364 FindThetaAtQuartz(thetaCerMax);
365 if(GetThetaAtQuartz() == 999.)
371 SetThetaPhotonInDRS(GetThetaAtQuartz());
372 SetPhiPhotonInDRS(phiPoint);
374 radiusMax = FromEmissionToCathode();
378 thetaCerMean = (thetaCerMax + thetaCerMin)/2;
380 FindThetaAtQuartz(thetaCerMean);
381 if(GetThetaAtQuartz() == 999.)
387 SetThetaPhotonInDRS(GetThetaAtQuartz());
388 SetPhiPhotonInDRS(phiPoint);
390 radiusMean = FromEmissionToCathode();
393 // AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax));
395 while (TMath::Abs(radiusMean-distPoint) > kTollerance)
398 if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean;
399 if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) {
401 thetaCerMin = thetaCerMean;
403 FindThetaAtQuartz(thetaCerMin);
404 SetThetaPhotonInDRS(GetThetaAtQuartz());
405 SetPhiPhotonInDRS(phiPoint);
407 radiusMin =FromEmissionToCathode();
410 thetaCerMean = (thetaCerMax + thetaCerMin)/2;
412 FindThetaAtQuartz(thetaCerMean);
413 SetThetaPhotonInDRS(GetThetaAtQuartz());
414 SetPhiPhotonInDRS(phiPoint);
416 radiusMean = FromEmissionToCathode();
420 // AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration));
421 SetThetaPhotonCerenkov(999.);
426 // AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean));
427 SetThetaPhotonCerenkov(thetaCerMean);
430 //__________________________________________________________________________________________________
431 void AliRICHRecon::FindAreaAndPortionOfRing()
433 //find fraction of the ring accepted by the RICH
435 Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing];
437 // Float_t xtoentr = GetEntranceX();
438 // Float_t ytoentr = GetEntranceY();
439 Float_t shiftX = GetShiftX();
440 Float_t shiftY = GetShiftY();
442 Float_t xemiss = GetXCoordOfEmission();
443 Float_t yemiss = GetYCoordOfEmission();
445 Float_t x0 = xemiss + shiftX;
446 Float_t y0 = yemiss + shiftY;
449 SetPhotonEnergy(AliRICHParam::MeanCkovEnergy());
450 SetFreonRefractiveIndex();
452 SetEmissionPoint(fRadiatorWidth/2.);
454 Float_t theta = GetThetaOfRing();
457 Int_t nPsiAccepted = 0;
460 for(Int_t i=0;i<NPointsOfRing-1;i++)
463 Float_t psi = 2*TMath::Pi()*i/NPointsOfRing;
465 SetThetaPhotonInTRS(theta);
466 SetPhiPhotonInTRS(psi);
467 FindPhotonAnglesInDRS();
469 Float_t radius = FromEmissionToCathode();
470 if (radius == 999.) continue;
474 Float_t xPointRing = GetXPointOnCathode() + shiftX;
475 Float_t yPointRing = GetYPointOnCathode() + shiftY;
477 SetDetectorWhereX(xPointRing);
478 SetDetectorWhereY(yPointRing);
480 Int_t zone = CheckDetectorAcceptance();
482 AliDebug(1,Form("acceptance to detector zone -> %d",zone));
486 FindIntersectionWithDetector();
487 xPoint[nPoints] = GetIntersectionX();
488 yPoint[nPoints] = GetIntersectionY();
492 xPoint[nPoints] = xPointRing;
493 yPoint[nPoints] = yPointRing;
501 xPoint[nPoints] = xPoint[0];
502 yPoint[nPoints] = yPoint[0];
508 for (Int_t i = 0; i < nPoints; i++)
510 area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0));
515 Float_t portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal));
519 SetPortionOfRing(portionOfRing);
521 //__________________________________________________________________________________________________
522 void AliRICHRecon::FindIntersectionWithDetector()
524 // find ring intersection with CsI edges
526 Float_t xIntersect, yIntersect;
527 Float_t x1, x2, y1, y2;
529 Float_t shiftX = GetShiftX();
530 Float_t shiftY = GetShiftY();
532 Float_t xPoint = GetXPointOnCathode() + shiftX;
533 Float_t yPoint = GetYPointOnCathode() + shiftY;
535 Float_t xemiss = GetXCoordOfEmission();
536 Float_t yemiss = GetYCoordOfEmission();
538 Float_t phi = GetPhiPhotonInDRS();
539 Float_t m = tan(phi);
541 Float_t x0 = xemiss + shiftX;
542 Float_t y0 = yemiss + shiftY;
566 yIntersect = m*(xIntersect - x0) + y0;
567 if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
569 SetIntersectionX(xIntersect);
570 SetIntersectionY(yIntersect);
575 yIntersect = m*(xIntersect - x0) + y0;
576 if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
578 SetIntersectionX(xIntersect);
579 SetIntersectionY(yIntersect);
584 xIntersect = (yIntersect - y0)/m + x0;
585 if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
587 SetIntersectionX(xIntersect);
588 SetIntersectionY(yIntersect);
593 xIntersect = (yIntersect - y0)/m + x0;
594 if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
596 SetIntersectionX(xIntersect);
597 SetIntersectionY(yIntersect);
601 cout << " sono fuori!!!!!!" << endl;
604 //__________________________________________________________________________________________________
605 Int_t AliRICHRecon::CheckDetectorAcceptance() const
607 // check for the acceptance
609 // crosses X -2.6 2.6 cm
612 Float_t xcoord = GetDetectorWhereX();
613 Float_t ycoord = GetDetectorWhereY();
617 if(ycoord > fYmax) return 2;
618 if(ycoord > fYmin && ycoord < fYmax) return 3;
619 if(ycoord < fYmin) return 4;
623 if(ycoord > fYmax) return 8;
624 if(ycoord > fYmin && ycoord < fYmax) return 7;
625 if(ycoord < fYmin) return 6;
627 if(xcoord > fXmin && xcoord < fXmax)
629 if(ycoord > fYmax) return 1;
630 if(ycoord > fYmin && ycoord < fYmax) return 0;
631 if(ycoord < fYmin) return 5;
635 //__________________________________________________________________________________________________
636 void AliRICHRecon::FindPhotonAnglesInDRS()
638 // Setup the rotation matrix of the track...
645 Float_t trackTheta = GetTrackTheta();
646 Float_t trackPhi = GetTrackPhi();
648 mtheta.RotateY(trackTheta);
649 mphi.RotateZ(trackPhi);
651 mrot = mphi * mtheta;
652 // minv = mrot.Inverse();
654 TVector3 photonInRadiator(1,1,1);
656 Float_t thetaCerenkov = GetThetaPhotonInTRS();
657 Float_t phiCerenkov = GetPhiPhotonInTRS();
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);
668 //__________________________________________________________________________________________________
669 Float_t AliRICHRecon::FromEmissionToCathode()
671 // trace from emission point to cathode
673 Float_t nfreon, nquartz, ngas;
675 SetFreonRefractiveIndex();
676 SetQuartzRefractiveIndex();
677 SetGasRefractiveIndex();
679 nfreon = GetFreonRefractiveIndex();
680 nquartz = GetQuartzRefractiveIndex();
681 ngas = GetGasRefractiveIndex();
683 Float_t trackTheta = GetTrackTheta();
684 Float_t trackPhi = GetTrackPhi();
685 Float_t lengthOfEmissionPoint = GetEmissionPoint();
687 Float_t theta = GetThetaPhotonInDRS();
688 Float_t phi = GetPhiPhotonInDRS();
690 // cout << " Theta " << Theta << " Phi " << Phi << endl;
692 Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi);
693 Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi);
695 SetXCoordOfEmission(xemiss);
696 SetYCoordOfEmission(yemiss);
698 Float_t thetaquar = SnellAngle(nfreon, nquartz, theta);
700 if(thetaquar == 999.)
702 SetXPointOnCathode(999.);
703 SetYPointOnCathode(999.);
707 Float_t thetagap = SnellAngle( nquartz, ngas, thetaquar);
711 SetXPointOnCathode(999.);
712 SetYPointOnCathode(999.);
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);
724 Float_t xtot = xemiss + xw + xq + xg;
725 Float_t ytot = yemiss + yw + yq + yg;
727 SetXPointOnCathode(xtot);
728 SetYPointOnCathode(ytot);
731 Float_t distanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2)
732 +TMath::Power(fPhotonLimitY,2));
734 return distanceFromEntrance;
737 //__________________________________________________________________________________________________
738 void AliRICHRecon::FindPhiPoint()
740 //find phi of generated point
742 Float_t xtoentr = GetEntranceX();
743 Float_t ytoentr = GetEntranceY();
745 Float_t trackTheta = GetTrackTheta();
746 Float_t trackPhi = GetTrackPhi();
748 Float_t emissionPoint = GetEmissionPoint();
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);
757 //__________________________________________________________________________________________________
758 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
760 // cerenkov angle from n and beta
762 // Compute the cerenkov angle
768 // cout << " warning in Cerenkoangle !!!!!! " << endl;
772 thetacer = acos (1./(n*beta));
775 //__________________________________________________________________________________________________
776 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
780 // Compute the Snell angle
782 Float_t sinrefractangle;
783 Float_t refractangle;
785 sinrefractangle = (n1/n2)*sin(theta1);
787 if(sinrefractangle>1.) {
788 // cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
793 refractangle = asin(sinrefractangle);
796 //__________________________________________________________________________________________________
797 void AliRICHRecon::HoughResponse()
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;
811 Int_t bin = (Int_t)(0.5+angle/(fDTheta));
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.;}
826 photsw->Fill(angle,weight);
827 SetPhotonWeight(weight);
829 for (Int_t i=1; i<=nBin;i++){
830 Int_t bin1= i-nCorrBand;
831 Int_t bin2= i+nCorrBand;
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);
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));
844 // Reset and delete objects
845 phots->Delete();photsw->Delete();resultw->Delete();
847 //__________________________________________________________________________________________________
848 void AliRICHRecon::FindThetaCerenkov()
850 // manage with weight for photons
853 Float_t weightThetaCerenkov = 0.;
855 Int_t nPhotons = GetPhotonsNumber();
856 for(Int_t i=0;i<nPhotons;i++)
860 if(GetPhotonFlag() == 2)
862 Float_t photonEta = GetPhotonEta();
863 Float_t photonWeight = GetPhotonWeight();
864 weightThetaCerenkov += photonEta*photonWeight;
871 weightThetaCerenkov /= wei;
875 weightThetaCerenkov = 0.;
878 SetThetaCerenkov(weightThetaCerenkov);
880 AliDebug(1,Form(" thetac weighted -> %f",weightThetaCerenkov));
882 //__________________________________________________________________________________________________
883 void AliRICHRecon::FlagPhotons()
887 Int_t nPhotonHough = 0;
889 Float_t thetaCerenkov = GetThetaCerenkov();
890 AliDebug(1,Form(" fThetaCerenkov %f ",thetaCerenkov));
892 Float_t thetaDist= thetaCerenkov - fThetaMin;
893 Int_t steps = (Int_t)(thetaDist / fDTheta);
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);
899 tmin = tavg - 0.5*fWindowWidth;
900 tmax = tavg + 0.5*fWindowWidth;
902 // Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
904 Int_t nPhotons = GetPhotonsNumber();
906 // for(Int_t i=0;i<candidatePhotonsNumber;i++)
908 for(Int_t i=0;i<nPhotons;i++)
912 Float_t photonEta = GetPhotonEta();
914 if(photonEta == -999.) continue;
916 if(photonEta >= tmin && photonEta <= tmax)
922 SetHoughPhotons(nPhotonHough);