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>
32 #include "AliRICHParam.h"
33 #include "AliRICHRecon.h"
34 #include "AliRICHHelix.h"
37 #define NPointsOfRing 201
39 //__________________________________________________________________________________________________
40 AliRICHRecon::AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId)
41 :TTask("RichRec","RichPat")
44 SetFreonScaleFactor(1);
46 if(pClusters->GetEntries()>200) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
48 fThetaMin = 0.0; fThetaMax = 0.75;
49 fDTheta = 0.001; fWindowWidth = 0.045;
51 fRadiatorWidth = AliRICHParam::Zfreon();
52 fQuartzWidth = AliRICHParam::Zwin();
53 fGapWidth = AliRICHParam::Freon2Pc() - fRadiatorWidth - fQuartzWidth;
55 fXmax = AliRICHParam::PcSizeX();
57 fYmax = AliRICHParam::PcSizeY();
58 SetTrackTheta(pHelix->Ploc().Theta());
59 SetTrackPhi(pHelix->Ploc().Phi());
61 SetShiftX(pHelix->PosRad().X());
62 SetShiftY(pHelix->PosRad().Y());
63 fpClusters = pClusters;
65 //__________________________________________________________________________________________________
66 Double_t AliRICHRecon::ThetaCerenkov()
68 // Pattern recognition method based on Hough transform
69 // Return theta Cerenkov for a given track and list of clusters which are set in ctor
71 if(fpClusters->GetEntries()==0) return -1;//no clusters at all for a given track
72 Bool_t kPatRec = kFALSE;
74 AliDebug(1,Form("---Track Parameters--- Theta: %f , Phi: %f ",GetTrackTheta()*TMath::RadToDeg(),GetTrackPhi()*TMath::RadToDeg()));
76 Int_t candidatePhotons = 0;
78 SetThetaCerenkov(999.);
80 SetHoughPhotonsNorm(0);
82 for (Int_t j=0; j < fpClusters->GetEntries(); j++){//clusters loop
87 if (j == GetMipIndex()) continue; // do not consider MIP cluster as a candidate photon
88 // if(((AliRICHCluster*)fpClusters->UncheckedAt(j))->Q()>AliRICHParam::QthMIP()) continue; // avoid MIP clusters from bkg
89 Float_t xtoentr = ((AliRICHCluster*)fpClusters->UncheckedAt(j))->X() - GetShiftX();
90 Float_t ytoentr = ((AliRICHCluster*)fpClusters->UncheckedAt(j))->Y() - GetShiftY();
91 SetEntranceX(xtoentr);
92 SetEntranceY(ytoentr);
94 // Int_t photonStatus = PhotonInBand();
95 // if(photonStatus == 0) continue;
97 FindThetaPhotonCerenkov();
98 Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
99 AliDebug(1,Form("THETA CERENKOV ---> %f",thetaPhotonCerenkov));
100 SetPhotonEta(thetaPhotonCerenkov);
104 if(candidatePhotons >= 1) kPatRec = kTRUE;
106 if(!kPatRec) return -1;
108 SetPhotonsNumber(fpClusters->GetEntries());
115 Int_t nPhotonHough = GetHoughPhotons();
119 SetThetaCerenkov(999.);
120 SetHoughPhotonsNorm(0.);
126 AliDebug(1,Form("Number of clusters accepted ---> %i",nPhotonHough));
128 // Float_t thetaCerenkov = GetThetaCerenkov();
129 // SetThetaOfRing(thetaCerenkov);
130 // FindAreaAndPortionOfRing();
132 // Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
133 // SetHoughPhotonsNorm(nPhotonHoughNorm);
135 // Calculate the area where the photon are accepted...
137 Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth;
138 SetThetaOfRing(thetaInternal);
139 FindAreaAndPortionOfRing();
140 Float_t internalArea = GetAreaOfRing();
142 Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth;
143 SetThetaOfRing(thetaExternal);
144 FindAreaAndPortionOfRing();
145 Float_t externalArea = GetAreaOfRing();
147 Float_t houghArea = externalArea - internalArea;
149 SetHoughArea(houghArea);
151 return GetThetaCerenkov();
154 //__________________________________________________________________________________________________
155 void AliRICHRecon::FindEmissionPoint()
157 //estimate the emission point in radiator
159 // Find emission point
161 Float_t absorbtionLenght=7.83*fRadiatorWidth; //absorption length in the freon (cm)
162 // 7.83 = -1/ln(T0) where
163 // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
164 Float_t photonLenght, photonLenghtMin, photonLenghtMax;
166 photonLenght=exp(-fRadiatorWidth/(absorbtionLenght*cos(fCerenkovAnglePad)));
167 photonLenghtMin=fRadiatorWidth*photonLenght/(1.-photonLenght);
168 photonLenghtMax=absorbtionLenght*cos(fCerenkovAnglePad);
169 Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax;
171 SetEmissionPoint(emissionPoint);
172 SetEmissionPoint(fRadiatorWidth/2); // tune the emission point
174 //__________________________________________________________________________________________________
175 Int_t AliRICHRecon::PhotonInBand()
177 //search band fro photon candidates
179 // Float_t massOfParticle;
184 Float_t xtoentr = GetEntranceX();
185 Float_t ytoentr = GetEntranceY();
190 Float_t phpad = GetPhiPoint();
194 SetPhotonEnergy(5.6);
195 SetEmissionPoint(fRadiatorWidth -0.0001);
196 SetFreonRefractiveIndex();
198 nfreon = GetFreonRefractiveIndex();
201 AliDebug(1,Form("thetacer in photoninband min %f",thetacer));
203 FindThetaAtQuartz(thetacer);
205 if(thetacer == 999. || GetThetaAtQuartz() == 999.)
208 SetXInnerRing(-999.);
209 SetYInnerRing(-999.);
210 SetRadiusInnerRing(-999.);
214 SetThetaPhotonInDRS(GetThetaAtQuartz());
215 SetPhiPhotonInDRS(phpad);
217 innerRadius = FromEmissionToCathode();
218 if(innerRadius == 999.) innerRadius = -999.;
220 SetXInnerRing(GetXPointOnCathode());
221 SetYInnerRing(GetYPointOnCathode());
222 SetRadiusInnerRing(innerRadius);
226 SetPhotonEnergy(7.7);
227 SetEmissionPoint(0.);
228 // SetMassHypotesis(0.139567);
229 SetFreonRefractiveIndex();
231 nfreon = GetFreonRefractiveIndex();
233 thetacer = Cerenkovangle(nfreon,1);
237 AliDebug(1,Form("thetacer in photoninband max %f",thetacer));
239 FindThetaAtQuartz(thetacer);
241 if(thetacer == 999. || GetThetaAtQuartz() == 999.)
246 SetRadiusOuterRing(999.);
250 SetThetaPhotonInDRS(GetThetaAtQuartz());
251 SetPhiPhotonInDRS(phpad);
253 outerRadius = FromEmissionToCathode();
254 // cout << " outerRadius " << outerRadius << endl;
255 SetXOuterRing(GetXPointOnCathode());
256 SetYOuterRing(GetYPointOnCathode());
257 SetRadiusOuterRing(outerRadius);
260 Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2));
262 AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius));
264 if(padradius>=innerRadius && padradius<=outerRadius) return 1;
267 //__________________________________________________________________________________________________
268 void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
270 //find the theta at the quartz plate
272 if(thetaCerenkov == 999.)
274 SetThetaAtQuartz(999.);
278 Float_t thetaAtQuartz = 999.;
280 Float_t trackTheta = GetTrackTheta();
282 if(trackTheta == 0) {
284 thetaAtQuartz = thetaCerenkov;
285 SetThetaAtQuartz(thetaAtQuartz);
289 Float_t trackPhi = GetTrackPhi();
290 Float_t phiPoint = GetPhiPoint();
292 Double_t den = TMath::Sin((Double_t)trackTheta)
293 *TMath::Cos((Double_t)trackPhi)
294 *TMath::Cos((Double_t)phiPoint) +
295 TMath::Sin((Double_t)trackTheta)
296 *TMath::Sin((Double_t)trackPhi)
297 *TMath::Sin((Double_t)phiPoint);
298 Double_t b = TMath::Cos((Double_t)trackTheta)/den;
299 Double_t c = -TMath::Cos((Double_t)thetaCerenkov)/den;
301 Double_t underSqrt = 1 + b*b - c*c;
304 SetThetaAtQuartz(999.);
308 Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c);
309 Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c);
311 Double_t thetaSol1 = 2*TMath::ATan(sol1);
312 Double_t thetaSol2 = 2*TMath::ATan(sol2);
314 if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1;
315 if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2;
317 // AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz));
319 SetThetaAtQuartz(thetaAtQuartz);
321 //__________________________________________________________________________________________________
322 void AliRICHRecon::FindThetaPhotonCerenkov()
324 //find theta cerenkov of ring
326 Float_t thetaCerMin = 0.;
327 Float_t thetaCerMax = 0.75;
328 Float_t thetaCerMean;
330 Float_t radiusMin, radiusMax, radiusMean;
331 Int_t nIteration = 0;
333 const Float_t kTollerance = 0.05;
336 Float_t phiPoint = GetPhiPoint();
338 SetPhotonEnergy(AliRICHParam::MeanCkovEnergy());
339 SetEmissionPoint(fRadiatorWidth/2);
341 Float_t xPoint = GetEntranceX();
342 Float_t yPoint = GetEntranceY();
343 Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint);
345 // AliDebug(1,Form(" DistPoint %f ",distPoint));
347 // Star minimization...
351 FindThetaAtQuartz(thetaCerMin);
353 if(GetThetaAtQuartz() == 999.)
359 SetThetaPhotonInDRS(GetThetaAtQuartz());
360 SetPhiPhotonInDRS(phiPoint);
362 radiusMin = FromEmissionToCathode();
367 FindThetaAtQuartz(thetaCerMax);
368 if(GetThetaAtQuartz() == 999.)
374 SetThetaPhotonInDRS(GetThetaAtQuartz());
375 SetPhiPhotonInDRS(phiPoint);
377 radiusMax = FromEmissionToCathode();
381 thetaCerMean = (thetaCerMax + thetaCerMin)/2;
383 FindThetaAtQuartz(thetaCerMean);
384 if(GetThetaAtQuartz() == 999.)
390 SetThetaPhotonInDRS(GetThetaAtQuartz());
391 SetPhiPhotonInDRS(phiPoint);
393 radiusMean = FromEmissionToCathode();
396 // AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax));
398 while (TMath::Abs(radiusMean-distPoint) > kTollerance)
401 if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean;
402 if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) {
404 thetaCerMin = thetaCerMean;
406 FindThetaAtQuartz(thetaCerMin);
407 SetThetaPhotonInDRS(GetThetaAtQuartz());
408 SetPhiPhotonInDRS(phiPoint);
410 radiusMin =FromEmissionToCathode();
413 thetaCerMean = (thetaCerMax + thetaCerMin)/2;
415 FindThetaAtQuartz(thetaCerMean);
416 SetThetaPhotonInDRS(GetThetaAtQuartz());
417 SetPhiPhotonInDRS(phiPoint);
419 radiusMean = FromEmissionToCathode();
423 // AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration));
424 SetThetaPhotonCerenkov(999.);
429 // AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean));
430 SetThetaPhotonCerenkov(thetaCerMean);
433 //__________________________________________________________________________________________________
434 void AliRICHRecon::FindAreaAndPortionOfRing()
436 //find fraction of the ring accepted by the RICH
438 Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing];
440 // Float_t xtoentr = GetEntranceX();
441 // Float_t ytoentr = GetEntranceY();
442 Float_t shiftX = GetShiftX();
443 Float_t shiftY = GetShiftY();
445 Float_t xemiss = GetXCoordOfEmission();
446 Float_t yemiss = GetYCoordOfEmission();
448 Float_t x0 = xemiss + shiftX;
449 Float_t y0 = yemiss + shiftY;
452 SetPhotonEnergy(AliRICHParam::MeanCkovEnergy());
453 SetFreonRefractiveIndex();
455 SetEmissionPoint(fRadiatorWidth/2.);
457 Float_t theta = GetThetaOfRing();
460 Int_t nPsiAccepted = 0;
463 for(Int_t i=0;i<NPointsOfRing-1;i++)
466 Float_t psi = 2*TMath::Pi()*i/NPointsOfRing;
468 SetThetaPhotonInTRS(theta);
469 SetPhiPhotonInTRS(psi);
470 FindPhotonAnglesInDRS();
472 Float_t radius = FromEmissionToCathode();
473 if (radius == 999.) continue;
477 Float_t xPointRing = GetXPointOnCathode() + shiftX;
478 Float_t yPointRing = GetYPointOnCathode() + shiftY;
480 SetDetectorWhereX(xPointRing);
481 SetDetectorWhereY(yPointRing);
483 Int_t zone = CheckDetectorAcceptance();
485 AliDebug(1,Form("acceptance to detector zone -> %d",zone));
489 FindIntersectionWithDetector();
490 xPoint[nPoints] = GetIntersectionX();
491 yPoint[nPoints] = GetIntersectionY();
495 xPoint[nPoints] = xPointRing;
496 yPoint[nPoints] = yPointRing;
504 xPoint[nPoints] = xPoint[0];
505 yPoint[nPoints] = yPoint[0];
511 for (Int_t i = 0; i < nPoints; i++)
513 area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0));
518 Float_t portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal));
522 SetPortionOfRing(portionOfRing);
524 //__________________________________________________________________________________________________
525 void AliRICHRecon::FindIntersectionWithDetector()
527 // find ring intersection with CsI edges
529 Float_t xIntersect, yIntersect;
530 Float_t x1, x2, y1, y2;
532 Float_t shiftX = GetShiftX();
533 Float_t shiftY = GetShiftY();
535 Float_t xPoint = GetXPointOnCathode() + shiftX;
536 Float_t yPoint = GetYPointOnCathode() + shiftY;
538 Float_t xemiss = GetXCoordOfEmission();
539 Float_t yemiss = GetYCoordOfEmission();
541 Float_t phi = GetPhiPhotonInDRS();
542 Float_t m = tan(phi);
544 Float_t x0 = xemiss + shiftX;
545 Float_t y0 = yemiss + shiftY;
569 yIntersect = m*(xIntersect - x0) + y0;
570 if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
572 SetIntersectionX(xIntersect);
573 SetIntersectionY(yIntersect);
578 yIntersect = m*(xIntersect - x0) + y0;
579 if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2)
581 SetIntersectionX(xIntersect);
582 SetIntersectionY(yIntersect);
587 xIntersect = (yIntersect - y0)/m + x0;
588 if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
590 SetIntersectionX(xIntersect);
591 SetIntersectionY(yIntersect);
596 xIntersect = (yIntersect - y0)/m + x0;
597 if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2)
599 SetIntersectionX(xIntersect);
600 SetIntersectionY(yIntersect);
604 cout << " sono fuori!!!!!!" << endl;
607 //__________________________________________________________________________________________________
608 Int_t AliRICHRecon::CheckDetectorAcceptance() const
610 // check for the acceptance
612 // crosses X -2.6 2.6 cm
615 Float_t xcoord = GetDetectorWhereX();
616 Float_t ycoord = GetDetectorWhereY();
620 if(ycoord > fYmax) return 2;
621 if(ycoord > fYmin && ycoord < fYmax) return 3;
622 if(ycoord < fYmin) return 4;
626 if(ycoord > fYmax) return 8;
627 if(ycoord > fYmin && ycoord < fYmax) return 7;
628 if(ycoord < fYmin) return 6;
630 if(xcoord > fXmin && xcoord < fXmax)
632 if(ycoord > fYmax) return 1;
633 if(ycoord > fYmin && ycoord < fYmax) return 0;
634 if(ycoord < fYmin) return 5;
638 //__________________________________________________________________________________________________
639 void AliRICHRecon::FindPhotonAnglesInDRS()
641 // Setup the rotation matrix of the track...
648 Float_t trackTheta = GetTrackTheta();
649 Float_t trackPhi = GetTrackPhi();
651 mtheta.RotateY(trackTheta);
652 mphi.RotateZ(trackPhi);
654 mrot = mphi * mtheta;
655 // minv = mrot.Inverse();
657 TVector3 photonInRadiator(1,1,1);
659 Float_t thetaCerenkov = GetThetaPhotonInTRS();
660 Float_t phiCerenkov = GetPhiPhotonInTRS();
662 photonInRadiator.SetTheta(thetaCerenkov);
663 photonInRadiator.SetPhi(phiCerenkov);
664 photonInRadiator = mrot * photonInRadiator;
665 Float_t theta = photonInRadiator.Theta();
666 Float_t phi = photonInRadiator.Phi();
667 SetThetaPhotonInDRS(theta);
668 SetPhiPhotonInDRS(phi);
671 //__________________________________________________________________________________________________
672 Float_t AliRICHRecon::FromEmissionToCathode()
674 // trace from emission point to cathode
676 Float_t nfreon, nquartz, ngas;
678 SetFreonRefractiveIndex();
679 SetQuartzRefractiveIndex();
680 SetGasRefractiveIndex();
682 nfreon = GetFreonRefractiveIndex();
683 nquartz = GetQuartzRefractiveIndex();
684 ngas = GetGasRefractiveIndex();
686 Float_t trackTheta = GetTrackTheta();
687 Float_t trackPhi = GetTrackPhi();
688 Float_t lengthOfEmissionPoint = GetEmissionPoint();
690 Float_t theta = GetThetaPhotonInDRS();
691 Float_t phi = GetPhiPhotonInDRS();
693 // cout << " Theta " << Theta << " Phi " << Phi << endl;
695 Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi);
696 Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi);
698 SetXCoordOfEmission(xemiss);
699 SetYCoordOfEmission(yemiss);
701 Float_t thetaquar = SnellAngle(nfreon, nquartz, theta);
703 if(thetaquar == 999.)
705 SetXPointOnCathode(999.);
706 SetYPointOnCathode(999.);
710 Float_t thetagap = SnellAngle( nquartz, ngas, thetaquar);
714 SetXPointOnCathode(999.);
715 SetYPointOnCathode(999.);
719 Float_t xw = (fRadiatorWidth - lengthOfEmissionPoint)*cos(phi)*tan(theta);
720 Float_t xq = fQuartzWidth*cos(phi)*tan(thetaquar);
721 Float_t xg = fGapWidth*cos(phi)*tan(thetagap);
722 Float_t yw = (fRadiatorWidth - lengthOfEmissionPoint)*sin(phi)*tan(theta);
723 Float_t yq = fQuartzWidth*sin(phi)*tan(thetaquar);
724 Float_t yg = fGapWidth*sin(phi)*tan(thetagap);
727 Float_t xtot = xemiss + xw + xq + xg;
728 Float_t ytot = yemiss + yw + yq + yg;
730 SetXPointOnCathode(xtot);
731 SetYPointOnCathode(ytot);
734 Float_t distanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2)
735 +TMath::Power(fPhotonLimitY,2));
737 return distanceFromEntrance;
740 //__________________________________________________________________________________________________
741 void AliRICHRecon::FindPhiPoint()
743 //find phi of generated point
745 Float_t xtoentr = GetEntranceX();
746 Float_t ytoentr = GetEntranceY();
748 Float_t trackTheta = GetTrackTheta();
749 Float_t trackPhi = GetTrackPhi();
751 Float_t emissionPoint = GetEmissionPoint();
753 Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi);
754 Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi);
755 Float_t phi = atan2(argY,argX);
760 //__________________________________________________________________________________________________
761 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
763 // cerenkov angle from n and beta
765 // Compute the cerenkov angle
771 // cout << " warning in Cerenkoangle !!!!!! " << endl;
775 thetacer = acos (1./(n*beta));
778 //__________________________________________________________________________________________________
779 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
783 // Compute the Snell angle
785 Float_t sinrefractangle;
786 Float_t refractangle;
788 sinrefractangle = (n1/n2)*sin(theta1);
790 if(sinrefractangle>1.) {
791 // cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
796 refractangle = asin(sinrefractangle);
799 //__________________________________________________________________________________________________
800 void AliRICHRecon::HoughResponse()
802 Int_t nChannels = (Int_t)(fThetaMax/fDTheta+0.5);
803 TH1F *phots = new TH1F("phots","phots",nChannels,0.,fThetaMax);
804 TH1F *photsw = new TH1F("photsw","photsw",nChannels,0.,fThetaMax);
805 TH1F *resultw = new TH1F("resultw","resultw",nChannels,0.,fThetaMax);
806 Int_t nPhotons = GetPhotonsNumber();
807 Int_t nBin = (Int_t)(fThetaMax/fDTheta);
808 Int_t nCorrBand = (Int_t)(fWindowWidth/(2*fDTheta));
809 AliDebug(1,Form("Ring reconstruction for track with theta %f",GetTrackTheta()*TMath::RadToDeg()));
810 for (Int_t kPhot=0; kPhot< nPhotons; kPhot++){
811 SetPhotonIndex(kPhot);
812 Double_t angle = GetPhotonEta();
813 if(angle<0||angle>fThetaMax) continue;
815 Int_t bin = (Int_t)(0.5+angle/(fDTheta));
818 Double_t lowerlimit = ((Float_t)bin)*fDTheta - 0.5*fDTheta;
819 SetThetaOfRing(lowerlimit);
820 FindAreaAndPortionOfRing();
821 Float_t area1 = GetAreaOfRing();
822 Double_t upperlimit = ((Float_t)bin)*fDTheta + 0.5*fDTheta;
823 SetThetaOfRing(upperlimit);
824 FindAreaAndPortionOfRing();
825 Float_t area2 = GetAreaOfRing();
826 AliDebug(1,Form("lowerlimit %f area %f ; upperlimit %f area %f",lowerlimit,area1,upperlimit,area2));
827 Float_t diffarea = area2 - area1;
828 if(diffarea>0){weight = 1./(area2-area1);}else{weight = 1.;}
830 AliDebug(1,Form("Calculated weight %f",weight));
831 photsw->Fill(angle,weight);
832 SetPhotonWeight(weight);
834 for (Int_t i=1; i<=nBin;i++){
835 Int_t bin1= i-nCorrBand;
836 Int_t bin2= i+nCorrBand;
838 if(bin2>nBin)bin2=nBin;
839 Double_t sumPhots=phots->Integral(bin1,bin2);
840 if(sumPhots<fMinNumPhots) continue; // cut on minimum n. of photons per ring
841 Double_t sumPhotsw=photsw->Integral(bin1,bin2);
842 resultw->Fill((Float_t)((i+0.5)*fDTheta),sumPhotsw);
844 // evaluate the "BEST" thetacerenkov....
845 Float_t *pVec = resultw->GetArray();
846 Int_t locMax = TMath::LocMax(nBin,pVec);
847 SetThetaCerenkov((Double_t)(locMax*fDTheta+0.5*fDTheta));
849 // Reset and delete objects
850 phots->Delete();photsw->Delete();resultw->Delete();
852 //__________________________________________________________________________________________________
853 void AliRICHRecon::FindThetaCerenkov()
855 // manage with weight for photons
858 Float_t weightThetaCerenkov = 0.;
860 Int_t nPhotons = GetPhotonsNumber();
861 Double_t etaMin=9999.,etaMax=0.;
862 for(Int_t i=0;i<nPhotons;i++)
866 if(GetPhotonFlag() == 2)
868 Float_t photonEta = GetPhotonEta();
869 if(photonEta<etaMin) etaMin=photonEta;
870 if(photonEta>etaMax) etaMax=photonEta;
871 Float_t photonWeight = GetPhotonWeight();
872 weightThetaCerenkov += photonEta*photonWeight;
877 if(wei != 0.) weightThetaCerenkov /= wei; else weightThetaCerenkov = 0.;
878 SetThetaCerenkov(weightThetaCerenkov);
880 // estimate of the n. of bkg photons
881 SetThetaOfRing(etaMin); FindAreaAndPortionOfRing(); Double_t internalArea = GetAreaOfRing();
882 SetThetaOfRing(etaMax); FindAreaAndPortionOfRing(); Double_t externalArea = GetAreaOfRing();
884 Double_t effArea = (AliRICHParam::PcSizeX()-AliRICHParam::DeadZone())*(AliRICHParam::PcSizeY()-2*AliRICHParam::DeadZone());
885 Double_t nPhotBKG = (externalArea-internalArea)/effArea*fpClusters->GetEntries();
886 if(nPhotBKG<0) nPhotBKG=0; //just protection from funny angles...
887 SetPhotBKG(nPhotBKG);
890 AliDebug(1,Form(" thetac weighted -> %f",weightThetaCerenkov));
892 //__________________________________________________________________________________________________
893 void AliRICHRecon::FlagPhotons()
897 Int_t nPhotonHough = 0;
899 Float_t thetaCerenkov = GetThetaCerenkov();
900 AliDebug(1,Form(" fThetaCerenkov %f ",thetaCerenkov));
902 Float_t thetaDist= thetaCerenkov - fThetaMin;
903 Int_t steps = (Int_t)(thetaDist / fDTheta);
905 Float_t tmin = fThetaMin + (Float_t)(steps - 1)*fDTheta;
906 Float_t tmax = fThetaMin + (Float_t)(steps)*fDTheta;
907 Float_t tavg = 0.5*(tmin+tmax);
909 tmin = tavg - 0.5*fWindowWidth;
910 tmax = tavg + 0.5*fWindowWidth;
912 // Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
914 Int_t nPhotons = GetPhotonsNumber();
916 // for(Int_t i=0;i<candidatePhotonsNumber;i++)
918 // TH1F* h1_phots = (TH1F*)gROOT->FindObjectAny("h1_phots");
920 // cout << "h1_phots " << h1_phots << endl;
922 for(Int_t i=0;i<nPhotons;i++)
926 Float_t photonEta = GetPhotonEta();
928 if(photonEta == -999.) continue;
930 if(photonEta >= tmin && photonEta <= tmax)
934 // if(h1_phots)h1_phots->Fill(photonEta);
937 SetHoughPhotons(nPhotonHough);