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")
42 fThetaMin = 0.0; fThetaMax = 0.75;
43 fDTheta = 0.001; fWindowWidth = 0.045;
45 fParam=AliRICHParam::Instance(); //get the pointer to AliRICHParam
47 //__________________________________________________________________________________________________
48 Double_t AliRICHRecon::ThetaCerenkov(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t &iMipId)
50 // Pattern recognition method based on Hough transform
51 // Return theta Cerenkov for a given track and list of clusters which are set in ctor
52 // 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.
54 // Returns: Track theta ckov in rad, nPhot contains number of photon candidates accepted for reconstruction track theta ckov
55 SetTrackTheta(pHelix->Ploc().Theta()); SetTrackPhi(pHelix->Ploc().Phi());
56 SetShiftX(pHelix->PosRad().X()); SetShiftY(pHelix->PosRad().Y());
57 fClusters = pClusters;
58 if(pClusters->GetEntries()>200) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
59 else fIsWEIGHT = kFALSE;
65 for (Int_t iClu=0; iClu<fClusters->GetEntriesFast();iClu++){//clusters loop
66 if(iClu == iMipId) continue; // do not consider MIP cluster as a photon candidate
71 AliRICHCluster *pClu=(AliRICHCluster*)fClusters->UncheckedAt(iClu); //get pointer to current cluster
72 // if(pClu->Q()>AliRICHParam::QthMIP()) continue; //avoid MIP clusters from bkg
73 SetEntranceX(pClu->X() - GetShiftX()); SetEntranceY(pClu->Y() - GetShiftY()); //cluster position with respect to track intersection
75 // if(PhotonInBand()==0) continue; ????????????
77 FindThetaPhotonCerenkov();
78 Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
79 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()
80 ,iClu,thetaPhotonCerenkov ));
81 SetPhotonEta(thetaPhotonCerenkov);
84 SetPhotonsNumber(fClusters->GetEntries());
86 if((iMipId=FlagPhotons(HoughResponse()))<1) return -11; //flag photons according to individual theta ckov with respect to most probable track theta ckov
89 // Float_t thetaCerenkov = GetThetaCerenkov();
90 // SetThetaOfRing(thetaCerenkov);
91 // FindAreaAndPortionOfRing();
93 // Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
94 // SetHoughPhotonsNorm(nPhotonHoughNorm);
96 // Calculate the area where the photon are accepted...
98 Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth;
99 SetThetaOfRing(thetaInternal);
100 FindAreaAndPortionOfRing();
101 Float_t internalArea = GetAreaOfRing();
103 Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth;
104 SetThetaOfRing(thetaExternal);
105 FindAreaAndPortionOfRing();
106 Float_t externalArea = GetAreaOfRing();
108 Float_t houghArea = externalArea - internalArea;
110 SetHoughArea(houghArea);
113 return GetThetaCerenkov();
115 //__________________________________________________________________________________________________
116 Int_t AliRICHRecon::PhotonInBand()
118 // Define valid band for photon candidates. For that photons with ThetaMin and ThetaMax are traced up to photcathode
124 Float_t xtoentr = GetEntranceX();
125 Float_t ytoentr = GetEntranceY();
130 Float_t phpad = GetPhiPoint();
134 SetEmissionPoint(AliRICHParam::RadThick() -0.0001);
136 nfreon = fParam->IdxC6F14(fParam->EckovMin());
139 AliDebug(1,Form("thetacer in photoninband min %f",thetacer));
141 FindThetaAtQuartz(thetacer);
143 if(thetacer == 999. || GetThetaAtQuartz() == 999.) {
145 SetXInnerRing(-999.);
146 SetYInnerRing(-999.);
147 SetRadiusInnerRing(-999.);
151 SetThetaPhotonInDRS(GetThetaAtQuartz());
152 SetPhiPhotonInDRS(phpad);
154 innerRadius = FromEmissionToCathode();
155 if(innerRadius == 999.) innerRadius = -999.;
157 SetXInnerRing(GetXPointOnCathode());
158 SetYInnerRing(GetYPointOnCathode());
159 SetRadiusInnerRing(innerRadius);
163 SetEmissionPoint(0.);
164 // SetMassHypotesis(0.139567);
166 nfreon = fParam->IdxC6F14(fParam->EckovMax());
168 thetacer = Cerenkovangle(nfreon,1);
172 AliDebug(1,Form("thetacer in photoninband max %f",thetacer));
174 FindThetaAtQuartz(thetacer);
176 if(thetacer == 999. || GetThetaAtQuartz() == 999.)
181 SetRadiusOuterRing(999.);
185 SetThetaPhotonInDRS(GetThetaAtQuartz());
186 SetPhiPhotonInDRS(phpad);
188 outerRadius = FromEmissionToCathode();
189 // cout << " outerRadius " << outerRadius << endl;
190 SetXOuterRing(GetXPointOnCathode());
191 SetYOuterRing(GetYPointOnCathode());
192 SetRadiusOuterRing(outerRadius);
195 Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2));
197 AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius));
199 if(padradius>=innerRadius && padradius<=outerRadius) return 1;
202 //__________________________________________________________________________________________________
203 void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
205 // find the theta at the quartz plate
207 if(thetaCerenkov == 999.) { SetThetaAtQuartz(999.); return; }
209 Float_t thetaAtQuartz = 999.;
211 Float_t trackTheta = GetTrackTheta();
213 if(trackTheta == 0) {
214 thetaAtQuartz = thetaCerenkov;
215 SetThetaAtQuartz(thetaAtQuartz);
219 Float_t trackPhi = GetTrackPhi();
220 Float_t phiPoint = GetPhiPoint();
222 Double_t den = TMath::Sin((Double_t)trackTheta)*TMath::Cos((Double_t)trackPhi)*TMath::Cos((Double_t)phiPoint) +
223 TMath::Sin((Double_t)trackTheta)*TMath::Sin((Double_t)trackPhi)*TMath::Sin((Double_t)phiPoint);
224 Double_t b = TMath::Cos((Double_t)trackTheta)/den;
225 Double_t c = -TMath::Cos((Double_t)thetaCerenkov)/den;
227 Double_t underSqrt = 1 + b*b - c*c;
230 SetThetaAtQuartz(999.);
234 Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c);
235 Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c);
237 Double_t thetaSol1 = 2*TMath::ATan(sol1);
238 Double_t thetaSol2 = 2*TMath::ATan(sol2);
240 if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1;
241 if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2;
243 // AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz));
245 SetThetaAtQuartz(thetaAtQuartz);
247 //__________________________________________________________________________________________________
248 void AliRICHRecon::FindThetaPhotonCerenkov()
250 //find theta cerenkov of ring
252 Float_t thetaCerMin = 0.;
253 Float_t thetaCerMax = 0.75;
254 Float_t thetaCerMean;
256 Float_t radiusMin, radiusMax, radiusMean;
257 Int_t nIteration = 0;
259 const Float_t kTollerance = 0.05;
262 Float_t phiPoint = GetPhiPoint();
264 SetEmissionPoint(AliRICHParam::RadThick()/2);
266 Float_t xPoint = GetEntranceX();
267 Float_t yPoint = GetEntranceY();
268 Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint);
270 // AliDebug(1,Form(" DistPoint %f ",distPoint));
272 // Star minimization...
276 FindThetaAtQuartz(thetaCerMin);
278 if(GetThetaAtQuartz() == 999.)
284 SetThetaPhotonInDRS(GetThetaAtQuartz());
285 SetPhiPhotonInDRS(phiPoint);
287 radiusMin = FromEmissionToCathode();
292 FindThetaAtQuartz(thetaCerMax);
293 if(GetThetaAtQuartz() == 999.)
299 SetThetaPhotonInDRS(GetThetaAtQuartz());
300 SetPhiPhotonInDRS(phiPoint);
302 radiusMax = FromEmissionToCathode();
306 thetaCerMean = (thetaCerMax + thetaCerMin)/2;
308 FindThetaAtQuartz(thetaCerMean);
309 if(GetThetaAtQuartz() == 999.)
315 SetThetaPhotonInDRS(GetThetaAtQuartz());
316 SetPhiPhotonInDRS(phiPoint);
318 radiusMean = FromEmissionToCathode();
321 // AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax));
323 while (TMath::Abs(radiusMean-distPoint) > kTollerance)
326 if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean;
327 if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) {
329 thetaCerMin = thetaCerMean;
331 FindThetaAtQuartz(thetaCerMin);
332 SetThetaPhotonInDRS(GetThetaAtQuartz());
333 SetPhiPhotonInDRS(phiPoint);
335 radiusMin =FromEmissionToCathode();
338 thetaCerMean = (thetaCerMax + thetaCerMin)/2;
340 FindThetaAtQuartz(thetaCerMean);
341 SetThetaPhotonInDRS(GetThetaAtQuartz());
342 SetPhiPhotonInDRS(phiPoint);
344 radiusMean = FromEmissionToCathode();
348 // AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration));
349 SetThetaPhotonCerenkov(999.);
354 // AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean));
355 SetThetaPhotonCerenkov(thetaCerMean);
358 //__________________________________________________________________________________________________
359 void AliRICHRecon::FindAreaAndPortionOfRing()
361 //find fraction of the ring accepted by the RICH
363 Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing];
365 // Float_t xtoentr = GetEntranceX();
366 // Float_t ytoentr = GetEntranceY();
367 Float_t shiftX = GetShiftX();
368 Float_t shiftY = GetShiftY();
370 Float_t xemiss = GetXCoordOfEmission();
371 Float_t yemiss = GetYCoordOfEmission();
373 Float_t x0 = xemiss + shiftX;
374 Float_t y0 = yemiss + shiftY;
378 SetEmissionPoint(AliRICHParam::RadThick()/2.);
380 Float_t theta = GetThetaOfRing();
383 Int_t nPsiAccepted = 0;
386 for(Int_t i=0;i<NPointsOfRing-1;i++){
387 Float_t psi = 2*TMath::Pi()*i/NPointsOfRing;
389 SetThetaPhotonInTRS(theta);
390 SetPhiPhotonInTRS(psi);
391 FindPhotonAnglesInDRS();
393 Float_t radius = FromEmissionToCathode();
394 if (radius == 999.) continue;
398 Float_t xPointRing = GetXPointOnCathode() + shiftX;
399 Float_t yPointRing = GetYPointOnCathode() + shiftY;
401 SetDetectorWhereX(xPointRing);
402 SetDetectorWhereY(yPointRing);
404 Int_t zone = CheckDetectorAcceptance();
406 AliDebug(1,Form("acceptance to detector zone -> %d",zone));
409 FindIntersectionWithDetector();
410 xPoint[nPoints] = GetIntersectionX(); yPoint[nPoints] = GetIntersectionY();
412 xPoint[nPoints] = xPointRing; yPoint[nPoints] = yPointRing;
418 xPoint[nPoints] = xPoint[0]; yPoint[nPoints] = yPoint[0];
424 for (Int_t i = 0; i < nPoints; i++)
426 area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0));
431 Float_t portionOfRing = 0;
433 portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal));
437 SetPortionOfRing(portionOfRing);
438 }//FindAreaAndPortionOfRing()
439 //__________________________________________________________________________________________________
440 void AliRICHRecon::FindIntersectionWithDetector()
442 // find ring intersection with CsI edges
444 Float_t xIntersect, yIntersect;
445 Float_t x1, x2, y1, y2;
447 Float_t shiftX = GetShiftX();
448 Float_t shiftY = GetShiftY();
450 Float_t xPoint = GetXPointOnCathode() + shiftX;
451 Float_t yPoint = GetYPointOnCathode() + shiftY;
453 Float_t xemiss = GetXCoordOfEmission();
454 Float_t yemiss = GetYCoordOfEmission();
456 Float_t phi = GetPhiPhotonInDRS();
457 Float_t m = tan(phi);
459 Float_t x0 = xemiss + shiftX;
460 Float_t y0 = yemiss + shiftY;
483 xIntersect = AliRICHParam::PcSizeX();
484 yIntersect = m*(xIntersect - x0) + y0;
485 if (yIntersect >= 0 && yIntersect <= AliRICHParam::PcSizeY() && xIntersect >= x1 && xIntersect <= x2)
487 SetIntersectionX(xIntersect);
488 SetIntersectionY(yIntersect);
493 yIntersect = m*(xIntersect - x0) + y0;
494 if (yIntersect >= 0 && yIntersect <= AliRICHParam::PcSizeY() && xIntersect >= x1 && xIntersect <= x2)
496 SetIntersectionX(xIntersect);
497 SetIntersectionY(yIntersect);
501 yIntersect = AliRICHParam::PcSizeY();
502 xIntersect = (yIntersect - y0)/m + x0;
503 if (xIntersect >= 0 && xIntersect <= AliRICHParam::PcSizeX() && yIntersect >= y1 && yIntersect <= y2)
505 SetIntersectionX(xIntersect);
506 SetIntersectionY(yIntersect);
511 xIntersect = (yIntersect - y0)/m + x0;
512 if (xIntersect >= 0 && xIntersect <= AliRICHParam::PcSizeX() && yIntersect >= y1 && yIntersect <= y2)
514 SetIntersectionX(xIntersect);
515 SetIntersectionY(yIntersect);
520 //__________________________________________________________________________________________________
521 Int_t AliRICHRecon::CheckDetectorAcceptance() const
523 // check for the acceptance
525 // crosses X -2.6 2.6 cm
528 Float_t xcoord = GetDetectorWhereX();
529 Float_t ycoord = GetDetectorWhereY();
531 if(xcoord > AliRICHParam::PcSizeX())
533 if(ycoord > AliRICHParam::PcSizeY()) return 2;
534 if(ycoord > 0 && ycoord < AliRICHParam::PcSizeY()) return 3;
535 if(ycoord < 0) return 4;
539 if(ycoord > AliRICHParam::PcSizeY()) return 8;
540 if(ycoord > 0 && ycoord < AliRICHParam::PcSizeY()) return 7;
541 if(ycoord < 0) return 6;
543 if(xcoord > 0 && xcoord < AliRICHParam::PcSizeX())
545 if(ycoord > AliRICHParam::PcSizeY()) return 1;
546 if(ycoord > 0 && ycoord < AliRICHParam::PcSizeY()) return 0;
547 if(ycoord < 0) return 5;
551 //__________________________________________________________________________________________________
552 void AliRICHRecon::FindPhotonAnglesInDRS()
554 // Setup the rotation matrix of the track...
561 Float_t trackTheta = GetTrackTheta();
562 Float_t trackPhi = GetTrackPhi();
564 mtheta.RotateY(trackTheta);
565 mphi.RotateZ(trackPhi);
567 mrot = mphi * mtheta;
568 // minv = mrot.Inverse();
570 TVector3 photonInRadiator(1,1,1);
572 Float_t thetaCerenkov = GetThetaPhotonInTRS();
573 Float_t phiCerenkov = GetPhiPhotonInTRS();
575 photonInRadiator.SetTheta(thetaCerenkov);
576 photonInRadiator.SetPhi(phiCerenkov);
577 photonInRadiator = mrot * photonInRadiator;
578 Float_t theta = photonInRadiator.Theta();
579 Float_t phi = photonInRadiator.Phi();
580 SetThetaPhotonInDRS(theta);
581 SetPhiPhotonInDRS(phi);
584 //__________________________________________________________________________________________________
585 Float_t AliRICHRecon::FromEmissionToCathode()
587 // Trace current photon from emission point somewhere in radiator to photocathode
591 Float_t nfreon, nquartz, ngas;
594 nfreon = fParam->IdxC6F14(fParam->EckovMean());
595 nquartz = fParam->IdxSiO2(fParam->EckovMean());
596 ngas = fParam->IdxCH4(fParam->EckovMean());
598 Float_t trackTheta = GetTrackTheta();
599 Float_t trackPhi = GetTrackPhi();
600 Float_t lengthOfEmissionPoint = GetEmissionPoint();
602 Float_t theta = GetThetaPhotonInDRS();
603 Float_t phi = GetPhiPhotonInDRS();
605 Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi);
606 Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi);
608 SetXCoordOfEmission(xemiss);
609 SetYCoordOfEmission(yemiss);
611 Float_t thetaquar = SnellAngle(nfreon, nquartz, theta);
613 if(thetaquar == 999.)
615 SetXPointOnCathode(999.);
616 SetYPointOnCathode(999.);
620 Float_t thetagap = SnellAngle( nquartz, ngas, thetaquar);
624 SetXPointOnCathode(999.);
625 SetYPointOnCathode(999.);
629 Float_t xw = (AliRICHParam::RadThick() - lengthOfEmissionPoint)*cos(phi)*tan(theta);
630 Float_t xq = AliRICHParam::WinThick()*cos(phi)*tan(thetaquar);
631 Float_t xg = AliRICHParam::Pc2Win()*cos(phi)*tan(thetagap);
632 Float_t yw = (AliRICHParam::RadThick() - lengthOfEmissionPoint)*sin(phi)*tan(theta);
633 Float_t yq = AliRICHParam::WinThick()*sin(phi)*tan(thetaquar);
634 Float_t yg = AliRICHParam::Pc2Win()*sin(phi)*tan(thetagap);
637 Float_t xtot = xemiss + xw + xq + xg;
638 Float_t ytot = yemiss + yw + yq + yg;
640 SetXPointOnCathode(xtot);
641 SetYPointOnCathode(ytot);
644 Float_t distanceFromEntrance = TMath::Sqrt(TMath::Power(fPhotonLimitX,2)+TMath::Power(fPhotonLimitY,2));
646 return distanceFromEntrance;
649 //__________________________________________________________________________________________________
650 void AliRICHRecon::FindPhiPoint()
652 //find phi of generated point
654 Float_t xtoentr = GetEntranceX();
655 Float_t ytoentr = GetEntranceY();
657 Float_t trackTheta = GetTrackTheta();
658 Float_t trackPhi = GetTrackPhi();
660 Float_t emissionPoint = GetEmissionPoint();
662 Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi);
663 Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi);
664 Float_t phi = atan2(argY,argX);
669 //__________________________________________________________________________________________________
670 Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta)
672 // cerenkov angle from n and beta
674 // Compute the cerenkov angle
680 // cout << " warning in Cerenkoangle !!!!!! " << endl;
684 thetacer = acos (1./(n*beta));
687 //__________________________________________________________________________________________________
688 Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
692 // Compute the Snell angle
694 Float_t sinrefractangle;
695 Float_t refractangle;
697 sinrefractangle = (n1/n2)*sin(theta1);
699 if(sinrefractangle>1.) {
700 // cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl;
705 refractangle = asin(sinrefractangle);
708 //__________________________________________________________________________________________________
709 Double_t AliRICHRecon::HoughResponse()
714 Int_t nChannels = (Int_t)(fThetaMax/fDTheta+0.5);
715 TH1F *phots = new TH1F("phots" ,"phots" ,nChannels,0,fThetaMax);
716 TH1F *photsw = new TH1F("photsw" ,"photsw" ,nChannels,0,fThetaMax);
717 TH1F *resultw = new TH1F("resultw","resultw",nChannels,0,fThetaMax);
718 Int_t nBin = (Int_t)(fThetaMax/fDTheta);
719 Int_t nCorrBand = (Int_t)(fWindowWidth/(2*fDTheta));
720 AliDebug(1,Form("Ring reconstruction for track with theta %f",GetTrackTheta()*TMath::RadToDeg()));
721 for (Int_t kPhot=0; kPhot< GetPhotonsNumber(); kPhot++){
722 SetPhotonIndex(kPhot);
723 Double_t angle = GetPhotonEta();
724 if(angle<0||angle>fThetaMax) continue;
726 Int_t bin = (Int_t)(0.5+angle/(fDTheta));
729 Double_t lowerlimit = ((Float_t)bin)*fDTheta - 0.5*fDTheta;
730 SetThetaOfRing(lowerlimit);
731 FindAreaAndPortionOfRing();
732 Float_t area1 = GetAreaOfRing();
733 Double_t upperlimit = ((Float_t)bin)*fDTheta + 0.5*fDTheta;
734 SetThetaOfRing(upperlimit);
735 FindAreaAndPortionOfRing();
736 Float_t area2 = GetAreaOfRing();
737 AliDebug(1,Form("lowerlimit %f area %f ; upperlimit %f area %f",lowerlimit,area1,upperlimit,area2));
738 Float_t diffarea = area2 - area1;
739 if(diffarea>0){weight = 1./(area2-area1);}else{weight = 1.;}
741 AliDebug(1,Form("Calculated weight %f",weight));
742 photsw->Fill(angle,weight);
743 SetPhotonWeight(weight);
745 for (Int_t i=1; i<=nBin;i++){
746 Int_t bin1= i-nCorrBand;
747 Int_t bin2= i+nCorrBand;
749 if(bin2>nBin)bin2=nBin;
750 Double_t sumPhots=phots->Integral(bin1,bin2);
751 if(sumPhots<fMinNumPhots) continue; // cut on minimum n. of photons per ring
752 Double_t sumPhotsw=photsw->Integral(bin1,bin2);
753 resultw->Fill((Float_t)((i+0.5)*fDTheta),sumPhotsw);
755 // evaluate the "BEST" theta ckov as the maximum value of histogramm
756 Float_t *pVec = resultw->GetArray();
757 Int_t locMax = TMath::LocMax(nBin,pVec);
758 phots->Delete();photsw->Delete();resultw->Delete(); // Reset and delete objects
760 return (Double_t)(locMax*fDTheta+0.5*fDTheta); //final most probable track theta ckov
762 //__________________________________________________________________________________________________
763 void AliRICHRecon::FindThetaCerenkov()
765 // manage with weight for photons
768 Float_t weightThetaCerenkov = 0.;
770 Double_t etaMin=9999.,etaMax=0.;
771 for(Int_t i=0;i<GetPhotonsNumber();i++){
773 if(GetPhotonFlag() == 2){
774 Float_t photonEta = GetPhotonEta();
775 if(photonEta<etaMin) etaMin=photonEta;
776 if(photonEta>etaMax) etaMax=photonEta;
777 Float_t photonWeight = GetPhotonWeight();
778 weightThetaCerenkov += photonEta*photonWeight;
783 if(wei != 0.) weightThetaCerenkov /= wei; else weightThetaCerenkov = 0.;
784 SetThetaCerenkov(weightThetaCerenkov);
786 // estimate of the n. of bkg photons
787 SetThetaOfRing(etaMin); FindAreaAndPortionOfRing(); Double_t internalArea = GetAreaOfRing();
788 SetThetaOfRing(etaMax); FindAreaAndPortionOfRing(); Double_t externalArea = GetAreaOfRing();
790 Double_t effArea = (AliRICHParam::PcSizeX()-AliRICHParam::DeadZone())*(AliRICHParam::PcSizeY()-2*AliRICHParam::DeadZone());
791 Double_t nPhotBKG = (externalArea-internalArea)/effArea*fClusters->GetEntries();
792 if(nPhotBKG<0) nPhotBKG=0; //just protection from funny angles...
793 SetPhotBKG(nPhotBKG);
796 AliDebug(1,Form(" thetac weighted -> %f",weightThetaCerenkov));
798 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
799 Int_t AliRICHRecon::FlagPhotons(Double_t thetaCkovHough)
801 // flag photon candidates if their individual theta ckov inside the window around theta ckov of Hough transform
802 // Arguments: thetaCkovHough- value of most probable theta ckov for track as returned by HoughResponse()
803 // Returns: number of photon candidates happened to be inside the window
805 Int_t steps = (Int_t)((thetaCkovHough - fThetaMin)/ fDTheta); //how many times we need to have fDTheta to fill the distance betwee fThetaMin and thetaCkovHough
807 Float_t tmin = fThetaMin + (Float_t)(steps - 1)*fDTheta;
808 Float_t tmax = fThetaMin + (Float_t)(steps)*fDTheta;
809 Float_t tavg = 0.5*(tmin+tmax);
811 tmin = tavg - 0.5*fWindowWidth; tmax = tavg + 0.5*fWindowWidth;
813 Int_t iInsideCnt = 0; //count photons which theta inside prdefined window
814 for(Int_t i=0;i<GetPhotonsNumber();i++){//photon candidates loop
815 SetPhotonIndex(i); Float_t photonEta = GetPhotonEta();
816 if(photonEta == -999.) continue;
817 if(photonEta >= tmin && photonEta <= tmax) { SetPhotonFlag(2); iInsideCnt++;}