]>
Commit | Line | Data |
---|---|---|
5cb4dfc3 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
88dd9ad4 | 16 | ////////////////////////////////////////////////////////////////////////// |
17 | // // | |
18 | // AliRICHRecon // | |
19 | // // | |
20 | // RICH class to perfom pattern recognition based on Hough transfrom // | |
21 | // // | |
22 | ////////////////////////////////////////////////////////////////////////// | |
23 | ||
5cb4dfc3 | 24 | #include <Riostream.h> |
5cb4dfc3 | 25 | #include <TMath.h> |
5cb4dfc3 | 26 | #include <TRotation.h> |
5cb4dfc3 | 27 | #include <TVector3.h> |
910735ae | 28 | #include <TH1F.h> |
024a7e64 | 29 | |
024a7e64 | 30 | #include "AliRICH.h" |
024a7e64 | 31 | #include "AliRICHParam.h" |
32 | #include "AliRICHRecon.h" | |
998b831f | 33 | #include "AliRICHHelix.h" |
34 | #include <AliLog.h> | |
5cb4dfc3 | 35 | |
36 | #define NPointsOfRing 201 | |
37 | ||
b068561d | 38 | //__________________________________________________________________________________________________ |
998b831f | 39 | AliRICHRecon::AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId) |
40 | :TTask("RichRec","RichPat") | |
5cb4dfc3 | 41 | { |
998b831f | 42 | // main ctor |
43 | SetFreonScaleFactor(1); | |
910735ae | 44 | // fIsWEIGHT = kFALSE; |
45 | fIsWEIGHT = kTRUE; | |
b068561d | 46 | fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75; |
e65514bb | 47 | fDTheta = 0.001; fWindowWidth = 0.045; |
910735ae | 48 | fMinNumPhots = 3; |
998b831f | 49 | fRadiatorWidth = AliRICHParam::Zfreon(); |
50 | fQuartzWidth = AliRICHParam::Zwin(); | |
51 | fGapWidth = AliRICHParam::Freon2Pc() - fRadiatorWidth - fQuartzWidth; | |
910735ae | 52 | fXmin = 0.; |
53 | fXmax = AliRICHParam::PcSizeX(); | |
54 | fYmin = 0.; | |
55 | fYmax = AliRICHParam::PcSizeY(); | |
998b831f | 56 | SetTrackTheta(pHelix->Ploc().Theta()); |
57 | SetTrackPhi(pHelix->Ploc().Phi()); | |
58 | SetMipIndex(iMipId); | |
59 | SetShiftX(pHelix->PosRad().X()); | |
60 | SetShiftY(pHelix->PosRad().Y()); | |
61 | fpClusters = pClusters; | |
5cb4dfc3 | 62 | } |
b068561d | 63 | //__________________________________________________________________________________________________ |
998b831f | 64 | Double_t AliRICHRecon::ThetaCerenkov() |
5cb4dfc3 | 65 | { |
998b831f | 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 | |
5cb4dfc3 | 68 | |
0fe8fa07 | 69 | if(fpClusters->GetEntries()==0) return -1;//no clusters at all for a given track |
5cb4dfc3 | 70 | Bool_t kPatRec = kFALSE; |
1cb5a97c | 71 | |
72 | AliDebug(1,Form("---Track Parameters--- Theta: %f , Phi: %f ",GetTrackTheta()*TMath::RadToDeg(),GetTrackPhi()*TMath::RadToDeg())); | |
5cb4dfc3 | 73 | |
88dd9ad4 | 74 | Int_t candidatePhotons = 0; |
5cb4dfc3 | 75 | |
5cb4dfc3 | 76 | SetThetaCerenkov(999.); |
77 | SetHoughPhotons(0); | |
78 | SetHoughPhotonsNorm(0); | |
5cb4dfc3 | 79 | |
998b831f | 80 | for (Int_t j=0; j < fpClusters->GetEntries(); j++){//clusters loop |
81 | SetPhotonIndex(j); | |
82 | SetPhotonFlag(0); | |
83 | SetPhotonEta(-999.); | |
84 | SetPhotonWeight(0.); | |
85 | if (j == GetMipIndex()) continue; // do not consider MIP cluster as a candidate photon | |
0fe8fa07 | 86 | Float_t xtoentr = ((AliRICHCluster*)fpClusters->UncheckedAt(j))->X() - GetShiftX(); |
87 | Float_t ytoentr = ((AliRICHCluster*)fpClusters->UncheckedAt(j))->Y() - GetShiftY(); | |
998b831f | 88 | SetEntranceX(xtoentr); |
89 | SetEntranceY(ytoentr); | |
90 | FindPhiPoint(); | |
91 | // Int_t photonStatus = PhotonInBand(); | |
92 | // if(photonStatus == 0) continue; | |
93 | SetPhotonFlag(1); | |
94 | FindThetaPhotonCerenkov(); | |
95 | Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov(); | |
1cb5a97c | 96 | AliDebug(1,Form("THETA CERENKOV ---> %f",thetaPhotonCerenkov)); |
998b831f | 97 | SetPhotonEta(thetaPhotonCerenkov); |
98 | candidatePhotons++; | |
99 | }//clusters loop | |
5cb4dfc3 | 100 | |
88dd9ad4 | 101 | if(candidatePhotons >= 1) kPatRec = kTRUE; |
5cb4dfc3 | 102 | |
0fe8fa07 | 103 | if(!kPatRec) return -1; |
998b831f | 104 | |
105 | SetPhotonsNumber(fpClusters->GetEntries()); | |
5cb4dfc3 | 106 | |
107 | HoughResponse(); | |
108 | ||
b068561d | 109 | fNrings++; |
5cb4dfc3 | 110 | |
111 | FlagPhotons(); | |
88dd9ad4 | 112 | Int_t nPhotonHough = GetHoughPhotons(); |
5cb4dfc3 | 113 | |
88dd9ad4 | 114 | if(nPhotonHough < 1) |
5cb4dfc3 | 115 | { |
116 | SetThetaCerenkov(999.); | |
117 | SetHoughPhotonsNorm(0.); | |
0fe8fa07 | 118 | return -1; |
5cb4dfc3 | 119 | } |
120 | ||
101624cd | 121 | FindThetaCerenkov(); |
5cb4dfc3 | 122 | |
998b831f | 123 | AliDebug(1,Form("Number of clusters accepted ---> %i",nPhotonHough)); |
124 | ||
101624cd | 125 | // Float_t thetaCerenkov = GetThetaCerenkov(); |
126 | // SetThetaOfRing(thetaCerenkov); | |
998b831f | 127 | // FindAreaAndPortionOfRing(); |
5cb4dfc3 | 128 | |
998b831f | 129 | // Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing(); |
130 | // SetHoughPhotonsNorm(nPhotonHoughNorm); | |
5cb4dfc3 | 131 | |
132 | // Calculate the area where the photon are accepted... | |
998b831f | 133 | /* |
88dd9ad4 | 134 | Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth; |
135 | SetThetaOfRing(thetaInternal); | |
5cb4dfc3 | 136 | FindAreaAndPortionOfRing(); |
88dd9ad4 | 137 | Float_t internalArea = GetAreaOfRing(); |
5cb4dfc3 | 138 | |
88dd9ad4 | 139 | Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth; |
140 | SetThetaOfRing(thetaExternal); | |
5cb4dfc3 | 141 | FindAreaAndPortionOfRing(); |
88dd9ad4 | 142 | Float_t externalArea = GetAreaOfRing(); |
5cb4dfc3 | 143 | |
88dd9ad4 | 144 | Float_t houghArea = externalArea - internalArea; |
5cb4dfc3 | 145 | |
88dd9ad4 | 146 | SetHoughArea(houghArea); |
998b831f | 147 | */ |
148 | return GetThetaCerenkov(); | |
5cb4dfc3 | 149 | |
998b831f | 150 | }//ThetaCerenkov() |
151 | //__________________________________________________________________________________________________ | |
5cb4dfc3 | 152 | void AliRICHRecon::FindEmissionPoint() |
88dd9ad4 | 153 | { |
154 | //estimate the emission point in radiator | |
5cb4dfc3 | 155 | |
156 | // Find emission point | |
157 | ||
88dd9ad4 | 158 | Float_t absorbtionLenght=7.83*fRadiatorWidth; //absorption length in the freon (cm) |
5cb4dfc3 | 159 | // 7.83 = -1/ln(T0) where |
160 | // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV) | |
161 | Float_t photonLenght, photonLenghtMin, photonLenghtMax; | |
162 | ||
88dd9ad4 | 163 | photonLenght=exp(-fRadiatorWidth/(absorbtionLenght*cos(fCerenkovAnglePad))); |
164 | photonLenghtMin=fRadiatorWidth*photonLenght/(1.-photonLenght); | |
5cb4dfc3 | 165 | photonLenghtMax=absorbtionLenght*cos(fCerenkovAnglePad); |
88dd9ad4 | 166 | Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax; |
5cb4dfc3 | 167 | |
88dd9ad4 | 168 | SetEmissionPoint(emissionPoint); |
998b831f | 169 | SetEmissionPoint(fRadiatorWidth/2); // tune the emission point |
5cb4dfc3 | 170 | } |
101624cd | 171 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 172 | Int_t AliRICHRecon::PhotonInBand() |
173 | { | |
88dd9ad4 | 174 | //search band fro photon candidates |
5cb4dfc3 | 175 | |
88dd9ad4 | 176 | // Float_t massOfParticle; |
5cb4dfc3 | 177 | Float_t nfreon; |
178 | ||
179 | Float_t thetacer; | |
180 | ||
88dd9ad4 | 181 | Float_t xtoentr = GetEntranceX(); |
182 | Float_t ytoentr = GetEntranceY(); | |
5cb4dfc3 | 183 | |
88dd9ad4 | 184 | Float_t innerRadius; |
185 | Float_t outerRadius; | |
5cb4dfc3 | 186 | |
187 | Float_t phpad = GetPhiPoint(); | |
188 | ||
5cb4dfc3 | 189 | |
190 | // inner radius // | |
191 | SetPhotonEnergy(5.6); | |
88dd9ad4 | 192 | SetEmissionPoint(fRadiatorWidth -0.0001); |
5cb4dfc3 | 193 | SetFreonRefractiveIndex(); |
194 | ||
5cb4dfc3 | 195 | nfreon = GetFreonRefractiveIndex(); |
5cb4dfc3 | 196 | thetacer = 0.; |
197 | ||
998b831f | 198 | AliDebug(1,Form("thetacer in photoninband min %f",thetacer)); |
5cb4dfc3 | 199 | |
200 | FindThetaAtQuartz(thetacer); | |
201 | ||
202 | if(thetacer == 999. || GetThetaAtQuartz() == 999.) | |
203 | { | |
88dd9ad4 | 204 | innerRadius = -999.; |
5cb4dfc3 | 205 | SetXInnerRing(-999.); |
206 | SetYInnerRing(-999.); | |
207 | SetRadiusInnerRing(-999.); | |
208 | } | |
209 | else | |
210 | { | |
211 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
212 | SetPhiPhotonInDRS(phpad); | |
213 | ||
88dd9ad4 | 214 | innerRadius = FromEmissionToCathode(); |
215 | if(innerRadius == 999.) innerRadius = -999.; | |
5cb4dfc3 | 216 | |
217 | SetXInnerRing(GetXPointOnCathode()); | |
218 | SetYInnerRing(GetYPointOnCathode()); | |
88dd9ad4 | 219 | SetRadiusInnerRing(innerRadius); |
5cb4dfc3 | 220 | } |
221 | ||
222 | // outer radius // | |
223 | SetPhotonEnergy(7.7); | |
224 | SetEmissionPoint(0.); | |
225 | // SetMassHypotesis(0.139567); | |
5cb4dfc3 | 226 | SetFreonRefractiveIndex(); |
227 | ||
5cb4dfc3 | 228 | nfreon = GetFreonRefractiveIndex(); |
229 | ||
998b831f | 230 | thetacer = Cerenkovangle(nfreon,1); |
5cb4dfc3 | 231 | |
232 | // thetacer = 0.75; | |
233 | ||
998b831f | 234 | AliDebug(1,Form("thetacer in photoninband max %f",thetacer)); |
5cb4dfc3 | 235 | |
236 | FindThetaAtQuartz(thetacer); | |
237 | ||
238 | if(thetacer == 999. || GetThetaAtQuartz() == 999.) | |
239 | { | |
88dd9ad4 | 240 | outerRadius = 999.; |
5cb4dfc3 | 241 | SetXOuterRing(999.); |
242 | SetYOuterRing(999.); | |
243 | SetRadiusOuterRing(999.); | |
244 | } | |
245 | else | |
246 | { | |
247 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
248 | SetPhiPhotonInDRS(phpad); | |
249 | ||
88dd9ad4 | 250 | outerRadius = FromEmissionToCathode(); |
251 | // cout << " outerRadius " << outerRadius << endl; | |
5cb4dfc3 | 252 | SetXOuterRing(GetXPointOnCathode()); |
253 | SetYOuterRing(GetYPointOnCathode()); | |
88dd9ad4 | 254 | SetRadiusOuterRing(outerRadius); |
5cb4dfc3 | 255 | } |
256 | ||
88dd9ad4 | 257 | Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2)); |
5cb4dfc3 | 258 | |
998b831f | 259 | AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius)); |
5cb4dfc3 | 260 | |
88dd9ad4 | 261 | if(padradius>=innerRadius && padradius<=outerRadius) return 1; |
5cb4dfc3 | 262 | return 0; |
263 | } | |
101624cd | 264 | //__________________________________________________________________________________________________ |
88dd9ad4 | 265 | void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov) |
5cb4dfc3 | 266 | { |
88dd9ad4 | 267 | //find the theta at the quartz plate |
5cb4dfc3 | 268 | |
88dd9ad4 | 269 | if(thetaCerenkov == 999.) |
5cb4dfc3 | 270 | { |
271 | SetThetaAtQuartz(999.); | |
272 | return; | |
273 | } | |
274 | ||
88dd9ad4 | 275 | Float_t thetaAtQuartz = 999.; |
5cb4dfc3 | 276 | |
88dd9ad4 | 277 | Float_t trackTheta = GetTrackTheta(); |
5cb4dfc3 | 278 | |
88dd9ad4 | 279 | if(trackTheta == 0) { |
5cb4dfc3 | 280 | |
88dd9ad4 | 281 | thetaAtQuartz = thetaCerenkov; |
282 | SetThetaAtQuartz(thetaAtQuartz); | |
5cb4dfc3 | 283 | return; |
284 | } | |
285 | ||
88dd9ad4 | 286 | Float_t trackPhi = GetTrackPhi(); |
287 | Float_t phiPoint = GetPhiPoint(); | |
5cb4dfc3 | 288 | |
88dd9ad4 | 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; | |
5cb4dfc3 | 297 | |
88dd9ad4 | 298 | Double_t underSqrt = 1 + b*b - c*c; |
5cb4dfc3 | 299 | |
88dd9ad4 | 300 | if(underSqrt < 0) { |
5cb4dfc3 | 301 | SetThetaAtQuartz(999.); |
302 | return; | |
303 | } | |
304 | ||
88dd9ad4 | 305 | Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c); |
306 | Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c); | |
5cb4dfc3 | 307 | |
88dd9ad4 | 308 | Double_t thetaSol1 = 2*TMath::ATan(sol1); |
309 | Double_t thetaSol2 = 2*TMath::ATan(sol2); | |
5cb4dfc3 | 310 | |
88dd9ad4 | 311 | if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1; |
312 | if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2; | |
5cb4dfc3 | 313 | |
998b831f | 314 | // AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz)); |
315 | ||
88dd9ad4 | 316 | SetThetaAtQuartz(thetaAtQuartz); |
5cb4dfc3 | 317 | } |
101624cd | 318 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 319 | void AliRICHRecon::FindThetaPhotonCerenkov() |
320 | { | |
88dd9ad4 | 321 | //find theta cerenkov of ring |
5cb4dfc3 | 322 | |
88dd9ad4 | 323 | Float_t thetaCerMin = 0.; |
324 | Float_t thetaCerMax = 0.75; | |
325 | Float_t thetaCerMean; | |
5cb4dfc3 | 326 | |
88dd9ad4 | 327 | Float_t radiusMin, radiusMax, radiusMean; |
5cb4dfc3 | 328 | Int_t nIteration = 0; |
329 | ||
88dd9ad4 | 330 | const Float_t kTollerance = 0.05; |
5cb4dfc3 | 331 | |
5cb4dfc3 | 332 | |
88dd9ad4 | 333 | Float_t phiPoint = GetPhiPoint(); |
5cb4dfc3 | 334 | |
101624cd | 335 | SetPhotonEnergy(AliRICHParam::MeanCkovEnergy()); |
88dd9ad4 | 336 | SetEmissionPoint(fRadiatorWidth/2); |
5cb4dfc3 | 337 | |
88dd9ad4 | 338 | Float_t xPoint = GetEntranceX(); |
339 | Float_t yPoint = GetEntranceY(); | |
998b831f | 340 | Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint); |
5cb4dfc3 | 341 | |
101624cd | 342 | // AliDebug(1,Form(" DistPoint %f ",distPoint)); |
5cb4dfc3 | 343 | |
344 | // Star minimization... | |
345 | ||
346 | // First value... | |
347 | ||
88dd9ad4 | 348 | FindThetaAtQuartz(thetaCerMin); |
5cb4dfc3 | 349 | |
350 | if(GetThetaAtQuartz() == 999.) | |
351 | { | |
88dd9ad4 | 352 | radiusMin = -999.; |
5cb4dfc3 | 353 | } |
354 | else | |
355 | { | |
356 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
88dd9ad4 | 357 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 358 | |
88dd9ad4 | 359 | radiusMin = FromEmissionToCathode(); |
5cb4dfc3 | 360 | } |
361 | ||
362 | // Second value... | |
363 | ||
88dd9ad4 | 364 | FindThetaAtQuartz(thetaCerMax); |
5cb4dfc3 | 365 | if(GetThetaAtQuartz() == 999.) |
366 | { | |
88dd9ad4 | 367 | radiusMax = 999.; |
5cb4dfc3 | 368 | } |
369 | else | |
370 | { | |
371 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
88dd9ad4 | 372 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 373 | |
88dd9ad4 | 374 | radiusMax = FromEmissionToCathode(); |
5cb4dfc3 | 375 | } |
376 | // Mean value... | |
377 | ||
88dd9ad4 | 378 | thetaCerMean = (thetaCerMax + thetaCerMin)/2; |
5cb4dfc3 | 379 | |
88dd9ad4 | 380 | FindThetaAtQuartz(thetaCerMean); |
5cb4dfc3 | 381 | if(GetThetaAtQuartz() == 999.) |
382 | { | |
88dd9ad4 | 383 | radiusMean = 999.; |
5cb4dfc3 | 384 | } |
385 | else | |
386 | { | |
387 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
88dd9ad4 | 388 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 389 | |
88dd9ad4 | 390 | radiusMean = FromEmissionToCathode(); |
5cb4dfc3 | 391 | } |
392 | ||
101624cd | 393 | // AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax)); |
5cb4dfc3 | 394 | |
88dd9ad4 | 395 | while (TMath::Abs(radiusMean-distPoint) > kTollerance) |
5cb4dfc3 | 396 | { |
397 | ||
88dd9ad4 | 398 | if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean; |
399 | if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) { | |
5cb4dfc3 | 400 | |
88dd9ad4 | 401 | thetaCerMin = thetaCerMean; |
5cb4dfc3 | 402 | |
88dd9ad4 | 403 | FindThetaAtQuartz(thetaCerMin); |
5cb4dfc3 | 404 | SetThetaPhotonInDRS(GetThetaAtQuartz()); |
88dd9ad4 | 405 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 406 | |
88dd9ad4 | 407 | radiusMin =FromEmissionToCathode(); |
5cb4dfc3 | 408 | } |
409 | ||
88dd9ad4 | 410 | thetaCerMean = (thetaCerMax + thetaCerMin)/2; |
5cb4dfc3 | 411 | |
88dd9ad4 | 412 | FindThetaAtQuartz(thetaCerMean); |
5cb4dfc3 | 413 | SetThetaPhotonInDRS(GetThetaAtQuartz()); |
88dd9ad4 | 414 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 415 | |
88dd9ad4 | 416 | radiusMean = FromEmissionToCathode(); |
5cb4dfc3 | 417 | |
418 | nIteration++; | |
419 | if(nIteration>=50) { | |
d0831219 | 420 | // AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration)); |
5cb4dfc3 | 421 | SetThetaPhotonCerenkov(999.); |
422 | return; | |
423 | } | |
424 | } | |
425 | ||
101624cd | 426 | // AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean)); |
88dd9ad4 | 427 | SetThetaPhotonCerenkov(thetaCerMean); |
5cb4dfc3 | 428 | |
429 | } | |
101624cd | 430 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 431 | void AliRICHRecon::FindAreaAndPortionOfRing() |
432 | { | |
88dd9ad4 | 433 | //find fraction of the ring accepted by the RICH |
5cb4dfc3 | 434 | |
88dd9ad4 | 435 | Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing]; |
5cb4dfc3 | 436 | |
88dd9ad4 | 437 | // Float_t xtoentr = GetEntranceX(); |
438 | // Float_t ytoentr = GetEntranceY(); | |
439 | Float_t shiftX = GetShiftX(); | |
440 | Float_t shiftY = GetShiftY(); | |
5cb4dfc3 | 441 | |
88dd9ad4 | 442 | Float_t xemiss = GetXCoordOfEmission(); |
443 | Float_t yemiss = GetYCoordOfEmission(); | |
5cb4dfc3 | 444 | |
88dd9ad4 | 445 | Float_t x0 = xemiss + shiftX; |
446 | Float_t y0 = yemiss + shiftY; | |
5cb4dfc3 | 447 | |
5cb4dfc3 | 448 | |
101624cd | 449 | SetPhotonEnergy(AliRICHParam::MeanCkovEnergy()); |
5cb4dfc3 | 450 | SetFreonRefractiveIndex(); |
451 | ||
88dd9ad4 | 452 | SetEmissionPoint(fRadiatorWidth/2.); |
5cb4dfc3 | 453 | |
88dd9ad4 | 454 | Float_t theta = GetThetaOfRing(); |
5cb4dfc3 | 455 | |
456 | Int_t nPoints = 0; | |
88dd9ad4 | 457 | Int_t nPsiAccepted = 0; |
458 | Int_t nPsiTotal = 0; | |
5cb4dfc3 | 459 | |
460 | for(Int_t i=0;i<NPointsOfRing-1;i++) | |
461 | { | |
462 | ||
88dd9ad4 | 463 | Float_t psi = 2*TMath::Pi()*i/NPointsOfRing; |
5cb4dfc3 | 464 | |
88dd9ad4 | 465 | SetThetaPhotonInTRS(theta); |
466 | SetPhiPhotonInTRS(psi); | |
5cb4dfc3 | 467 | FindPhotonAnglesInDRS(); |
468 | ||
88dd9ad4 | 469 | Float_t radius = FromEmissionToCathode(); |
470 | if (radius == 999.) continue; | |
5cb4dfc3 | 471 | |
88dd9ad4 | 472 | nPsiTotal++; |
5cb4dfc3 | 473 | |
88dd9ad4 | 474 | Float_t xPointRing = GetXPointOnCathode() + shiftX; |
475 | Float_t yPointRing = GetYPointOnCathode() + shiftY; | |
5cb4dfc3 | 476 | |
88dd9ad4 | 477 | SetDetectorWhereX(xPointRing); |
478 | SetDetectorWhereY(yPointRing); | |
5cb4dfc3 | 479 | |
88dd9ad4 | 480 | Int_t zone = CheckDetectorAcceptance(); |
5cb4dfc3 | 481 | |
910735ae | 482 | AliDebug(1,Form("acceptance to detector zone -> %d",zone)); |
5cb4dfc3 | 483 | |
88dd9ad4 | 484 | if (zone != 0) |
5cb4dfc3 | 485 | { |
486 | FindIntersectionWithDetector(); | |
88dd9ad4 | 487 | xPoint[nPoints] = GetIntersectionX(); |
488 | yPoint[nPoints] = GetIntersectionY(); | |
5cb4dfc3 | 489 | } |
490 | else | |
491 | { | |
88dd9ad4 | 492 | xPoint[nPoints] = xPointRing; |
493 | yPoint[nPoints] = yPointRing; | |
494 | nPsiAccepted++; | |
5cb4dfc3 | 495 | } |
496 | ||
497 | nPoints++; | |
498 | ||
499 | } | |
500 | ||
88dd9ad4 | 501 | xPoint[nPoints] = xPoint[0]; |
502 | yPoint[nPoints] = yPoint[0]; | |
5cb4dfc3 | 503 | |
504 | // find area... | |
505 | ||
88dd9ad4 | 506 | Float_t area = 0; |
5cb4dfc3 | 507 | |
508 | for (Int_t i = 0; i < nPoints; i++) | |
509 | { | |
88dd9ad4 | 510 | area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0)); |
5cb4dfc3 | 511 | } |
512 | ||
88dd9ad4 | 513 | area *= 0.5; |
5cb4dfc3 | 514 | |
88dd9ad4 | 515 | Float_t portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal)); |
5cb4dfc3 | 516 | |
5cb4dfc3 | 517 | |
88dd9ad4 | 518 | SetAreaOfRing(area); |
519 | SetPortionOfRing(portionOfRing); | |
5cb4dfc3 | 520 | } |
101624cd | 521 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 522 | void AliRICHRecon::FindIntersectionWithDetector() |
523 | { | |
88dd9ad4 | 524 | // find ring intersection with CsI edges |
5cb4dfc3 | 525 | |
88dd9ad4 | 526 | Float_t xIntersect, yIntersect; |
5cb4dfc3 | 527 | Float_t x1, x2, y1, y2; |
528 | ||
88dd9ad4 | 529 | Float_t shiftX = GetShiftX(); |
530 | Float_t shiftY = GetShiftY(); | |
5cb4dfc3 | 531 | |
88dd9ad4 | 532 | Float_t xPoint = GetXPointOnCathode() + shiftX; |
533 | Float_t yPoint = GetYPointOnCathode() + shiftY; | |
5cb4dfc3 | 534 | |
88dd9ad4 | 535 | Float_t xemiss = GetXCoordOfEmission(); |
536 | Float_t yemiss = GetYCoordOfEmission(); | |
5cb4dfc3 | 537 | |
88dd9ad4 | 538 | Float_t phi = GetPhiPhotonInDRS(); |
539 | Float_t m = tan(phi); | |
5cb4dfc3 | 540 | |
88dd9ad4 | 541 | Float_t x0 = xemiss + shiftX; |
542 | Float_t y0 = yemiss + shiftY; | |
5cb4dfc3 | 543 | |
88dd9ad4 | 544 | if(xPoint > x0) |
5cb4dfc3 | 545 | { |
546 | x1 = x0; | |
88dd9ad4 | 547 | x2 = xPoint; |
5cb4dfc3 | 548 | } |
549 | else | |
550 | { | |
551 | x2 = x0; | |
88dd9ad4 | 552 | x1 = xPoint; |
5cb4dfc3 | 553 | } |
88dd9ad4 | 554 | if(yPoint > y0) |
5cb4dfc3 | 555 | { |
556 | y1 = y0; | |
88dd9ad4 | 557 | y2 = yPoint; |
5cb4dfc3 | 558 | } |
559 | else | |
560 | { | |
561 | y2 = y0; | |
88dd9ad4 | 562 | y1 = yPoint; |
5cb4dfc3 | 563 | } |
564 | // | |
88dd9ad4 | 565 | xIntersect = fXmax; |
566 | yIntersect = m*(xIntersect - x0) + y0; | |
567 | if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2) | |
5cb4dfc3 | 568 | { |
88dd9ad4 | 569 | SetIntersectionX(xIntersect); |
570 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 571 | return; |
572 | } | |
573 | // | |
88dd9ad4 | 574 | xIntersect = fXmin; |
575 | yIntersect = m*(xIntersect - x0) + y0; | |
576 | if (yIntersect >= fYmin && yIntersect <= fYmax && xIntersect >= x1 && xIntersect <= x2) | |
5cb4dfc3 | 577 | { |
88dd9ad4 | 578 | SetIntersectionX(xIntersect); |
579 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 580 | return; |
581 | } | |
582 | // | |
88dd9ad4 | 583 | yIntersect = fYmax; |
584 | xIntersect = (yIntersect - y0)/m + x0; | |
585 | if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2) | |
5cb4dfc3 | 586 | { |
88dd9ad4 | 587 | SetIntersectionX(xIntersect); |
588 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 589 | return; |
590 | } | |
591 | // | |
88dd9ad4 | 592 | yIntersect = fYmin; |
593 | xIntersect = (yIntersect - y0)/m + x0; | |
594 | if (xIntersect >= fXmin && xIntersect <= fXmax && yIntersect >= y1 && yIntersect <= y2) | |
5cb4dfc3 | 595 | { |
88dd9ad4 | 596 | SetIntersectionX(xIntersect); |
597 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 598 | return; |
599 | } | |
600 | ||
601 | cout << " sono fuori!!!!!!" << endl; | |
5cb4dfc3 | 602 | |
603 | } | |
b068561d | 604 | //__________________________________________________________________________________________________ |
88dd9ad4 | 605 | Int_t AliRICHRecon::CheckDetectorAcceptance() const |
5cb4dfc3 | 606 | { |
88dd9ad4 | 607 | // check for the acceptance |
5cb4dfc3 | 608 | |
609 | // crosses X -2.6 2.6 cm | |
610 | // crosses Y -1 1 cm | |
611 | ||
88dd9ad4 | 612 | Float_t xcoord = GetDetectorWhereX(); |
613 | Float_t ycoord = GetDetectorWhereY(); | |
5cb4dfc3 | 614 | |
88dd9ad4 | 615 | if(xcoord > fXmax) |
5cb4dfc3 | 616 | { |
88dd9ad4 | 617 | if(ycoord > fYmax) return 2; |
618 | if(ycoord > fYmin && ycoord < fYmax) return 3; | |
619 | if(ycoord < fYmin) return 4; | |
5cb4dfc3 | 620 | } |
88dd9ad4 | 621 | if(xcoord < fXmin) |
5cb4dfc3 | 622 | { |
88dd9ad4 | 623 | if(ycoord > fYmax) return 8; |
624 | if(ycoord > fYmin && ycoord < fYmax) return 7; | |
625 | if(ycoord < fYmin) return 6; | |
5cb4dfc3 | 626 | } |
88dd9ad4 | 627 | if(xcoord > fXmin && xcoord < fXmax) |
5cb4dfc3 | 628 | { |
88dd9ad4 | 629 | if(ycoord > fYmax) return 1; |
630 | if(ycoord > fYmin && ycoord < fYmax) return 0; | |
631 | if(ycoord < fYmin) return 5; | |
5cb4dfc3 | 632 | } |
633 | return 999; | |
634 | } | |
b068561d | 635 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 636 | void AliRICHRecon::FindPhotonAnglesInDRS() |
637 | { | |
638 | // Setup the rotation matrix of the track... | |
639 | ||
88dd9ad4 | 640 | TRotation mtheta; |
641 | TRotation mphi; | |
642 | TRotation minv; | |
643 | TRotation mrot; | |
5cb4dfc3 | 644 | |
88dd9ad4 | 645 | Float_t trackTheta = GetTrackTheta(); |
646 | Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 647 | |
88dd9ad4 | 648 | mtheta.RotateY(trackTheta); |
649 | mphi.RotateZ(trackPhi); | |
5cb4dfc3 | 650 | |
88dd9ad4 | 651 | mrot = mphi * mtheta; |
652 | // minv = mrot.Inverse(); | |
5cb4dfc3 | 653 | |
88dd9ad4 | 654 | TVector3 photonInRadiator(1,1,1); |
5cb4dfc3 | 655 | |
88dd9ad4 | 656 | Float_t thetaCerenkov = GetThetaPhotonInTRS(); |
657 | Float_t phiCerenkov = GetPhiPhotonInTRS(); | |
5cb4dfc3 | 658 | |
88dd9ad4 | 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); | |
5cb4dfc3 | 666 | |
667 | } | |
101624cd | 668 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 669 | Float_t AliRICHRecon::FromEmissionToCathode() |
670 | { | |
88dd9ad4 | 671 | // trace from emission point to cathode |
5cb4dfc3 | 672 | |
673 | Float_t nfreon, nquartz, ngas; | |
674 | ||
675 | SetFreonRefractiveIndex(); | |
676 | SetQuartzRefractiveIndex(); | |
677 | SetGasRefractiveIndex(); | |
678 | ||
679 | nfreon = GetFreonRefractiveIndex(); | |
680 | nquartz = GetQuartzRefractiveIndex(); | |
681 | ngas = GetGasRefractiveIndex(); | |
682 | ||
88dd9ad4 | 683 | Float_t trackTheta = GetTrackTheta(); |
684 | Float_t trackPhi = GetTrackPhi(); | |
685 | Float_t lengthOfEmissionPoint = GetEmissionPoint(); | |
5cb4dfc3 | 686 | |
88dd9ad4 | 687 | Float_t theta = GetThetaPhotonInDRS(); |
688 | Float_t phi = GetPhiPhotonInDRS(); | |
5cb4dfc3 | 689 | |
690 | // cout << " Theta " << Theta << " Phi " << Phi << endl; | |
691 | ||
88dd9ad4 | 692 | Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi); |
693 | Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi); | |
5cb4dfc3 | 694 | |
88dd9ad4 | 695 | SetXCoordOfEmission(xemiss); |
696 | SetYCoordOfEmission(yemiss); | |
5cb4dfc3 | 697 | |
88dd9ad4 | 698 | Float_t thetaquar = SnellAngle(nfreon, nquartz, theta); |
5cb4dfc3 | 699 | |
700 | if(thetaquar == 999.) | |
701 | { | |
702 | SetXPointOnCathode(999.); | |
703 | SetYPointOnCathode(999.); | |
704 | return thetaquar; | |
705 | } | |
706 | ||
707 | Float_t thetagap = SnellAngle( nquartz, ngas, thetaquar); | |
708 | ||
709 | if(thetagap == 999.) | |
710 | { | |
711 | SetXPointOnCathode(999.); | |
712 | SetYPointOnCathode(999.); | |
713 | return thetagap; | |
714 | } | |
715 | ||
88dd9ad4 | 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); | |
5cb4dfc3 | 722 | |
5cb4dfc3 | 723 | |
88dd9ad4 | 724 | Float_t xtot = xemiss + xw + xq + xg; |
725 | Float_t ytot = yemiss + yw + yq + yg; | |
5cb4dfc3 | 726 | |
727 | SetXPointOnCathode(xtot); | |
728 | SetYPointOnCathode(ytot); | |
729 | ||
5cb4dfc3 | 730 | |
88dd9ad4 | 731 | Float_t distanceFromEntrance = sqrt(TMath::Power(fPhotonLimitX,2) |
5cb4dfc3 | 732 | +TMath::Power(fPhotonLimitY,2)); |
733 | ||
88dd9ad4 | 734 | return distanceFromEntrance; |
5cb4dfc3 | 735 | |
736 | } | |
101624cd | 737 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 738 | void AliRICHRecon::FindPhiPoint() |
739 | { | |
88dd9ad4 | 740 | //find phi of generated point |
5cb4dfc3 | 741 | |
88dd9ad4 | 742 | Float_t xtoentr = GetEntranceX(); |
743 | Float_t ytoentr = GetEntranceY(); | |
5cb4dfc3 | 744 | |
88dd9ad4 | 745 | Float_t trackTheta = GetTrackTheta(); |
746 | Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 747 | |
88dd9ad4 | 748 | Float_t emissionPoint = GetEmissionPoint(); |
5cb4dfc3 | 749 | |
88dd9ad4 | 750 | Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi); |
751 | Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi); | |
101624cd | 752 | Float_t phi = atan2(argY,argX); |
5cb4dfc3 | 753 | |
998b831f | 754 | SetPhiPoint(phi); |
5cb4dfc3 | 755 | |
756 | } | |
101624cd | 757 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 758 | Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta) |
759 | { | |
88dd9ad4 | 760 | // cerenkov angle from n and beta |
5cb4dfc3 | 761 | |
762 | // Compute the cerenkov angle | |
763 | ||
764 | Float_t thetacer; | |
765 | ||
766 | if((n*beta)<1.) { | |
767 | thetacer = 999.; | |
768 | // cout << " warning in Cerenkoangle !!!!!! " << endl; | |
769 | return thetacer; | |
770 | } | |
771 | ||
772 | thetacer = acos (1./(n*beta)); | |
773 | return thetacer; | |
774 | } | |
101624cd | 775 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 776 | Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1) |
88dd9ad4 | 777 | { |
778 | // Snell law | |
5cb4dfc3 | 779 | |
780 | // Compute the Snell angle | |
781 | ||
782 | Float_t sinrefractangle; | |
783 | Float_t refractangle; | |
784 | ||
785 | sinrefractangle = (n1/n2)*sin(theta1); | |
786 | ||
787 | if(sinrefractangle>1.) { | |
788 | // cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl; | |
789 | refractangle = 999.; | |
790 | return refractangle; | |
791 | } | |
792 | ||
793 | refractangle = asin(sinrefractangle); | |
794 | return refractangle; | |
795 | } | |
101624cd | 796 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 797 | void AliRICHRecon::HoughResponse() |
88dd9ad4 | 798 | { |
910735ae | 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); | |
88dd9ad4 | 802 | Int_t nPhotons = GetPhotonsNumber(); |
910735ae | 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(); | |
5cb4dfc3 | 809 | if(angle == -999.) continue; |
910735ae | 810 | phots->Fill(angle); |
811 | Int_t bin = (Int_t)(0.5+angle/(fDTheta)); | |
812 | Double_t weight=1.; | |
813 | if(fIsWEIGHT){ | |
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.;} | |
825 | } | |
826 | photsw->Fill(angle,weight); | |
827 | SetPhotonWeight(weight); | |
828 | } | |
829 | for (Int_t i=1; i<=nBin;i++){ | |
830 | Int_t bin1= i-nCorrBand; | |
831 | Int_t bin2= i+nCorrBand; | |
832 | if(bin1<1) bin1=1; | |
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); | |
838 | } | |
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)); | |
843 | ||
844 | // Reset and delete objects | |
845 | phots->Delete();photsw->Delete();resultw->Delete(); | |
846 | }//HoughResponse | |
101624cd | 847 | //__________________________________________________________________________________________________ |
848 | void AliRICHRecon::FindThetaCerenkov() | |
5cb4dfc3 | 849 | { |
88dd9ad4 | 850 | // manage with weight for photons |
5cb4dfc3 | 851 | |
852 | Float_t wei = 0.; | |
88dd9ad4 | 853 | Float_t weightThetaCerenkov = 0.; |
5cb4dfc3 | 854 | |
88dd9ad4 | 855 | Int_t nPhotons = GetPhotonsNumber(); |
856 | for(Int_t i=0;i<nPhotons;i++) | |
5cb4dfc3 | 857 | { |
858 | SetPhotonIndex(i); | |
859 | ||
860 | if(GetPhotonFlag() == 2) | |
861 | { | |
88dd9ad4 | 862 | Float_t photonEta = GetPhotonEta(); |
863 | Float_t photonWeight = GetPhotonWeight(); | |
864 | weightThetaCerenkov += photonEta*photonWeight; | |
865 | wei += photonWeight; | |
5cb4dfc3 | 866 | } |
867 | } | |
868 | ||
869 | if(wei != 0.) | |
870 | { | |
88dd9ad4 | 871 | weightThetaCerenkov /= wei; |
5cb4dfc3 | 872 | } |
873 | else | |
874 | { | |
88dd9ad4 | 875 | weightThetaCerenkov = 0.; |
5cb4dfc3 | 876 | } |
877 | ||
88dd9ad4 | 878 | SetThetaCerenkov(weightThetaCerenkov); |
5cb4dfc3 | 879 | |
101624cd | 880 | AliDebug(1,Form(" thetac weighted -> %f",weightThetaCerenkov)); |
5cb4dfc3 | 881 | } |
101624cd | 882 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 883 | void AliRICHRecon::FlagPhotons() |
884 | { | |
88dd9ad4 | 885 | // flag photons |
5cb4dfc3 | 886 | |
88dd9ad4 | 887 | Int_t nPhotonHough = 0; |
5cb4dfc3 | 888 | |
88dd9ad4 | 889 | Float_t thetaCerenkov = GetThetaCerenkov(); |
1cb5a97c | 890 | AliDebug(1,Form(" fThetaCerenkov %f ",thetaCerenkov)); |
5cb4dfc3 | 891 | |
88dd9ad4 | 892 | Float_t thetaDist= thetaCerenkov - fThetaMin; |
893 | Int_t steps = (Int_t)(thetaDist / fDTheta); | |
5cb4dfc3 | 894 | |
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); | |
898 | ||
899 | tmin = tavg - 0.5*fWindowWidth; | |
900 | tmax = tavg + 0.5*fWindowWidth; | |
901 | ||
88dd9ad4 | 902 | // Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber(); |
5cb4dfc3 | 903 | |
88dd9ad4 | 904 | Int_t nPhotons = GetPhotonsNumber(); |
5cb4dfc3 | 905 | |
88dd9ad4 | 906 | // for(Int_t i=0;i<candidatePhotonsNumber;i++) |
5cb4dfc3 | 907 | |
88dd9ad4 | 908 | for(Int_t i=0;i<nPhotons;i++) |
5cb4dfc3 | 909 | { |
910 | SetPhotonIndex(i); | |
911 | ||
88dd9ad4 | 912 | Float_t photonEta = GetPhotonEta(); |
5cb4dfc3 | 913 | |
88dd9ad4 | 914 | if(photonEta == -999.) continue; |
5cb4dfc3 | 915 | |
88dd9ad4 | 916 | if(photonEta >= tmin && photonEta <= tmax) |
5cb4dfc3 | 917 | { |
918 | SetPhotonFlag(2); | |
88dd9ad4 | 919 | nPhotonHough++; |
5cb4dfc3 | 920 | } |
921 | } | |
88dd9ad4 | 922 | SetHoughPhotons(nPhotonHough); |
5cb4dfc3 | 923 | } |