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