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