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