]>
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 // | |
db910db9 | 21 | // for single chamber // |
88dd9ad4 | 22 | ////////////////////////////////////////////////////////////////////////// |
23 | ||
db910db9 | 24 | #include "AliRICHRecon.h" //class header |
5cb4dfc3 | 25 | #include <TMath.h> |
5cb4dfc3 | 26 | #include <TRotation.h> |
5cb4dfc3 | 27 | #include <TVector3.h> |
910735ae | 28 | #include <TH1F.h> |
024a7e64 | 29 | |
db910db9 | 30 | #include "AliRICHCluster.h" //ThetaCerenkov() |
024a7e64 | 31 | #include "AliRICHParam.h" |
db910db9 | 32 | #include "AliRICHHelix.h" //ThetaCerenkov() |
998b831f | 33 | #include <AliLog.h> |
5cb4dfc3 | 34 | |
35 | #define NPointsOfRing 201 | |
36 | ||
b068561d | 37 | //__________________________________________________________________________________________________ |
5c09947e | 38 | AliRICHRecon::AliRICHRecon(): |
39 | TTask ("RichRec","RichPat"), | |
5c09947e | 40 | fPhotonsNumber(0), |
41 | fPhotonIndex(0), // should we use -1? | |
42 | fFittedHoughPhotons(0), | |
43 | fMinNumPhots(3), | |
44 | fTrackTheta(0), | |
45 | fTrackPhi(0), | |
46 | fMinDist(999), | |
47 | fTrackBeta(0), | |
48 | fXtoentr(999), | |
49 | fYtoentr(999), | |
50 | fThetaPhotonInTRS(0), | |
51 | fPhiPhotonInTRS(0), | |
52 | fThetaPhotonInDRS(0), | |
53 | fPhiPhotonInDRS(0), | |
54 | fThetaAtQuartz(0), | |
55 | fXEmiss(-999), | |
56 | fYEmiss(-999), | |
57 | fXInner(-999), | |
58 | fYInner(-999), | |
59 | fXOuter(-999), | |
60 | fYOuter(-999), | |
61 | fInnerRadius(-999), | |
62 | fOuterRadius(-999), | |
63 | fEphot(0), | |
64 | fLengthEmissionPoint(0), | |
65 | fPhotonLimitX(999), | |
66 | fPhotonLimitY(999), | |
67 | fDistanceFromCluster(999), | |
68 | fCerenkovAnglePad(999), | |
69 | fThetaPhotonCerenkov(0), | |
70 | fShiftX(0), | |
71 | fShiftY(0), | |
72 | fXcoord(999), | |
73 | fYcoord(999), | |
74 | fIntersectionX(999), | |
75 | fIntersectionY(999), | |
76 | fMassHypotesis(0.139567), | |
77 | fThetaOfRing(0), | |
78 | fAreaOfRing(0), | |
79 | fPortionOfRing(0), | |
80 | fHoughArea(0), | |
81 | fHoughRMS(999), | |
5c09947e | 82 | fFittedTrackTheta(0), |
83 | fFittedTrackPhi(0), | |
84 | fFittedThetaCerenkov(0), | |
85 | fThetaMin(0.0), | |
86 | fThetaMax(0.75), | |
87 | fIsWEIGHT(kFALSE), | |
88 | fIsBACKGROUND(kFALSE), | |
89 | fDTheta(0.001), | |
90 | fWindowWidth(0.045), | |
91 | fNumEtaPhotons(0), | |
92 | fnPhotBKG(0), | |
93 | fThetaCerenkov(0), | |
94 | fWeightThetaCerenkov(0), | |
4c3a1f92 | 95 | fThetaPeakPos(0), |
96 | fRingSigma2(0) | |
5cb4dfc3 | 97 | { |
998b831f | 98 | // main ctor |
5c09947e | 99 | for (Int_t i=0; i<3000; i++) { |
100 | fPhotonFlag[i] = 0; | |
101 | fPhiPoint[i] = -999; | |
102 | fPhotonEta[i] = -999; | |
103 | fPhotonWeight[i] = 0.0; | |
104 | fEtaFlag[i] = 0; | |
105 | fEtaPhotons[i] = 0; | |
106 | fWeightPhotons[i] = 0; | |
107 | } | |
5cb4dfc3 | 108 | } |
b068561d | 109 | //__________________________________________________________________________________________________ |
aad5d262 | 110 | Double_t AliRICHRecon::ThetaCerenkov(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t &nphot) |
5cb4dfc3 | 111 | { |
998b831f | 112 | // Pattern recognition method based on Hough transform |
113 | // Return theta Cerenkov for a given track and list of clusters which are set in ctor | |
db910db9 | 114 | // Remeber that list of clusters must contain more then 1 cluster. This considiration implies that normally we have 1 mip cluster and few photon clusters per track. |
115 | // Argume | |
aad5d262 | 116 | // Returns: Track theta ckov in rad, nphot is the number of photon candidates accepted for reconstruction track theta ckov |
db910db9 | 117 | SetTrackTheta(pHelix->Ploc().Theta()); SetTrackPhi(pHelix->Ploc().Phi()); |
118 | SetShiftX(pHelix->PosRad().X()); SetShiftY(pHelix->PosRad().Y()); | |
db910db9 | 119 | if(pClusters->GetEntries()>200) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction |
120 | else fIsWEIGHT = kFALSE; | |
5cb4dfc3 | 121 | |
5cb4dfc3 | 122 | |
5cb4dfc3 | 123 | |
db910db9 | 124 | SetThetaCerenkov(-1); |
5cb4dfc3 | 125 | |
84093f6f | 126 | // |
127 | // Photon Flag: Flag = 0 initial set; Flag = 1 good candidate (charge compatible with photon); Flag = 2 photon used for the ring; | |
128 | // | |
129 | ||
4c3a1f92 | 130 | for (Int_t iClu=0; iClu<pClusters->GetEntriesFast();iClu++){//clusters loop |
aad5d262 | 131 | if(iClu == nphot) continue; // do not consider MIP cluster as a photon candidate |
db910db9 | 132 | SetPhotonIndex(iClu); |
998b831f | 133 | SetPhotonFlag(0); |
134 | SetPhotonEta(-999.); | |
135 | SetPhotonWeight(0.); | |
4c3a1f92 | 136 | AliRICHCluster *pClu=(AliRICHCluster*)pClusters->UncheckedAt(iClu); //get pointer to current cluster |
84093f6f | 137 | if(pClu->Q()>AliRICHParam::QthMIP()) continue; //avoid MIP clusters from bkg |
db910db9 | 138 | SetEntranceX(pClu->X() - GetShiftX()); SetEntranceY(pClu->Y() - GetShiftY()); //cluster position with respect to track intersection |
998b831f | 139 | FindPhiPoint(); |
84093f6f | 140 | SetPhotonFlag(1); |
998b831f | 141 | FindThetaPhotonCerenkov(); |
142 | Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov(); | |
db910db9 | 143 | AliDebug(1,Form("Track Theta=%5.2f deg, Phi=%5.2f deg Photon clus=%2i ThetaCkov=%5.2f rad",GetTrackTheta()*TMath::RadToDeg(),GetTrackPhi()*TMath::RadToDeg() |
144 | ,iClu,thetaPhotonCerenkov )); | |
998b831f | 145 | SetPhotonEta(thetaPhotonCerenkov); |
998b831f | 146 | }//clusters loop |
5cb4dfc3 | 147 | |
4c3a1f92 | 148 | SetPhotonsNumber(pClusters->GetEntries()); |
998b831f | 149 | |
aad5d262 | 150 | if((nphot=FlagPhotons(HoughResponse()))<1) return -11; //flag photons according to individual theta ckov with respect to most probable track theta ckov |
5cb4dfc3 | 151 | |
5cb4dfc3 | 152 | |
101624cd | 153 | // Float_t thetaCerenkov = GetThetaCerenkov(); |
154 | // SetThetaOfRing(thetaCerenkov); | |
998b831f | 155 | // FindAreaAndPortionOfRing(); |
5cb4dfc3 | 156 | |
998b831f | 157 | // Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing(); |
158 | // SetHoughPhotonsNorm(nPhotonHoughNorm); | |
5cb4dfc3 | 159 | |
160 | // Calculate the area where the photon are accepted... | |
998b831f | 161 | /* |
88dd9ad4 | 162 | Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth; |
163 | SetThetaOfRing(thetaInternal); | |
5cb4dfc3 | 164 | FindAreaAndPortionOfRing(); |
88dd9ad4 | 165 | Float_t internalArea = GetAreaOfRing(); |
5cb4dfc3 | 166 | |
88dd9ad4 | 167 | Float_t thetaExternal = thetaCerenkov + 0.5*fWindowWidth; |
168 | SetThetaOfRing(thetaExternal); | |
5cb4dfc3 | 169 | FindAreaAndPortionOfRing(); |
88dd9ad4 | 170 | Float_t externalArea = GetAreaOfRing(); |
5cb4dfc3 | 171 | |
88dd9ad4 | 172 | Float_t houghArea = externalArea - internalArea; |
5cb4dfc3 | 173 | |
88dd9ad4 | 174 | SetHoughArea(houghArea); |
998b831f | 175 | */ |
4c3a1f92 | 176 | FindThetaCerenkov(pClusters->GetEntries()); |
998b831f | 177 | return GetThetaCerenkov(); |
998b831f | 178 | }//ThetaCerenkov() |
179 | //__________________________________________________________________________________________________ | |
5cb4dfc3 | 180 | Int_t AliRICHRecon::PhotonInBand() |
181 | { | |
db910db9 | 182 | // Define valid band for photon candidates. For that photons with ThetaMin and ThetaMax are traced up to photcathode |
5cb4dfc3 | 183 | |
5cb4dfc3 | 184 | Float_t nfreon; |
185 | ||
186 | Float_t thetacer; | |
187 | ||
88dd9ad4 | 188 | Float_t xtoentr = GetEntranceX(); |
189 | Float_t ytoentr = GetEntranceY(); | |
5cb4dfc3 | 190 | |
88dd9ad4 | 191 | Float_t innerRadius; |
192 | Float_t outerRadius; | |
5cb4dfc3 | 193 | |
194 | Float_t phpad = GetPhiPoint(); | |
195 | ||
5cb4dfc3 | 196 | |
197 | // inner radius // | |
db910db9 | 198 | SetEmissionPoint(AliRICHParam::RadThick() -0.0001); |
5cb4dfc3 | 199 | |
4c3a1f92 | 200 | nfreon = AliRICHParam::Instance()->IdxC6F14(AliRICHParam::Instance()->EckovMin()); |
5cb4dfc3 | 201 | thetacer = 0.; |
202 | ||
998b831f | 203 | AliDebug(1,Form("thetacer in photoninband min %f",thetacer)); |
5cb4dfc3 | 204 | |
205 | FindThetaAtQuartz(thetacer); | |
206 | ||
db910db9 | 207 | if(thetacer == 999. || GetThetaAtQuartz() == 999.) { |
88dd9ad4 | 208 | innerRadius = -999.; |
5cb4dfc3 | 209 | SetXInnerRing(-999.); |
210 | SetYInnerRing(-999.); | |
211 | SetRadiusInnerRing(-999.); | |
212 | } | |
213 | else | |
214 | { | |
215 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
216 | SetPhiPhotonInDRS(phpad); | |
217 | ||
88dd9ad4 | 218 | innerRadius = FromEmissionToCathode(); |
219 | if(innerRadius == 999.) innerRadius = -999.; | |
5cb4dfc3 | 220 | |
221 | SetXInnerRing(GetXPointOnCathode()); | |
222 | SetYInnerRing(GetYPointOnCathode()); | |
88dd9ad4 | 223 | SetRadiusInnerRing(innerRadius); |
5cb4dfc3 | 224 | } |
225 | ||
226 | // outer radius // | |
5cb4dfc3 | 227 | SetEmissionPoint(0.); |
228 | // SetMassHypotesis(0.139567); | |
5cb4dfc3 | 229 | |
4c3a1f92 | 230 | nfreon = AliRICHParam::Instance()->IdxC6F14(AliRICHParam::Instance()->EckovMax()); |
5cb4dfc3 | 231 | |
998b831f | 232 | thetacer = Cerenkovangle(nfreon,1); |
5cb4dfc3 | 233 | |
234 | // thetacer = 0.75; | |
235 | ||
998b831f | 236 | AliDebug(1,Form("thetacer in photoninband max %f",thetacer)); |
5cb4dfc3 | 237 | |
238 | FindThetaAtQuartz(thetacer); | |
239 | ||
240 | if(thetacer == 999. || GetThetaAtQuartz() == 999.) | |
241 | { | |
88dd9ad4 | 242 | outerRadius = 999.; |
5cb4dfc3 | 243 | SetXOuterRing(999.); |
244 | SetYOuterRing(999.); | |
245 | SetRadiusOuterRing(999.); | |
246 | } | |
247 | else | |
248 | { | |
249 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
250 | SetPhiPhotonInDRS(phpad); | |
251 | ||
88dd9ad4 | 252 | outerRadius = FromEmissionToCathode(); |
253 | // cout << " outerRadius " << outerRadius << endl; | |
5cb4dfc3 | 254 | SetXOuterRing(GetXPointOnCathode()); |
255 | SetYOuterRing(GetYPointOnCathode()); | |
88dd9ad4 | 256 | SetRadiusOuterRing(outerRadius); |
5cb4dfc3 | 257 | } |
258 | ||
88dd9ad4 | 259 | Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2)); |
5cb4dfc3 | 260 | |
998b831f | 261 | AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius)); |
5cb4dfc3 | 262 | |
88dd9ad4 | 263 | if(padradius>=innerRadius && padradius<=outerRadius) return 1; |
5cb4dfc3 | 264 | return 0; |
db910db9 | 265 | }//PhotonInBand() |
101624cd | 266 | //__________________________________________________________________________________________________ |
88dd9ad4 | 267 | void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov) |
5cb4dfc3 | 268 | { |
db910db9 | 269 | // find the theta at the quartz plate |
5cb4dfc3 | 270 | |
db910db9 | 271 | if(thetaCerenkov == 999.) { SetThetaAtQuartz(999.); return; } |
5cb4dfc3 | 272 | |
88dd9ad4 | 273 | Float_t thetaAtQuartz = 999.; |
5cb4dfc3 | 274 | |
88dd9ad4 | 275 | Float_t trackTheta = GetTrackTheta(); |
5cb4dfc3 | 276 | |
88dd9ad4 | 277 | if(trackTheta == 0) { |
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 | |
db910db9 | 286 | Double_t den = TMath::Sin((Double_t)trackTheta)*TMath::Cos((Double_t)trackPhi)*TMath::Cos((Double_t)phiPoint) + |
287 | TMath::Sin((Double_t)trackTheta)*TMath::Sin((Double_t)trackPhi)*TMath::Sin((Double_t)phiPoint); | |
88dd9ad4 | 288 | Double_t b = TMath::Cos((Double_t)trackTheta)/den; |
289 | Double_t c = -TMath::Cos((Double_t)thetaCerenkov)/den; | |
5cb4dfc3 | 290 | |
88dd9ad4 | 291 | Double_t underSqrt = 1 + b*b - c*c; |
5cb4dfc3 | 292 | |
88dd9ad4 | 293 | if(underSqrt < 0) { |
5cb4dfc3 | 294 | SetThetaAtQuartz(999.); |
295 | return; | |
296 | } | |
297 | ||
88dd9ad4 | 298 | Double_t sol1 = (1+TMath::Sqrt(underSqrt))/(b-c); |
299 | Double_t sol2 = (1-TMath::Sqrt(underSqrt))/(b-c); | |
5cb4dfc3 | 300 | |
88dd9ad4 | 301 | Double_t thetaSol1 = 2*TMath::ATan(sol1); |
302 | Double_t thetaSol2 = 2*TMath::ATan(sol2); | |
5cb4dfc3 | 303 | |
88dd9ad4 | 304 | if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1; |
305 | if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2; | |
5cb4dfc3 | 306 | |
998b831f | 307 | // AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz)); |
308 | ||
88dd9ad4 | 309 | SetThetaAtQuartz(thetaAtQuartz); |
5cb4dfc3 | 310 | } |
101624cd | 311 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 312 | void AliRICHRecon::FindThetaPhotonCerenkov() |
313 | { | |
88dd9ad4 | 314 | //find theta cerenkov of ring |
5cb4dfc3 | 315 | |
88dd9ad4 | 316 | Float_t thetaCerMin = 0.; |
317 | Float_t thetaCerMax = 0.75; | |
318 | Float_t thetaCerMean; | |
5cb4dfc3 | 319 | |
88dd9ad4 | 320 | Float_t radiusMin, radiusMax, radiusMean; |
5cb4dfc3 | 321 | Int_t nIteration = 0; |
322 | ||
88dd9ad4 | 323 | const Float_t kTollerance = 0.05; |
5cb4dfc3 | 324 | |
5cb4dfc3 | 325 | |
88dd9ad4 | 326 | Float_t phiPoint = GetPhiPoint(); |
5cb4dfc3 | 327 | |
db910db9 | 328 | SetEmissionPoint(AliRICHParam::RadThick()/2); |
5cb4dfc3 | 329 | |
88dd9ad4 | 330 | Float_t xPoint = GetEntranceX(); |
331 | Float_t yPoint = GetEntranceY(); | |
998b831f | 332 | Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint); |
5cb4dfc3 | 333 | |
101624cd | 334 | // AliDebug(1,Form(" DistPoint %f ",distPoint)); |
5cb4dfc3 | 335 | |
336 | // Star minimization... | |
337 | ||
338 | // First value... | |
339 | ||
88dd9ad4 | 340 | FindThetaAtQuartz(thetaCerMin); |
5cb4dfc3 | 341 | |
342 | if(GetThetaAtQuartz() == 999.) | |
343 | { | |
88dd9ad4 | 344 | radiusMin = -999.; |
5cb4dfc3 | 345 | } |
346 | else | |
347 | { | |
348 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
88dd9ad4 | 349 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 350 | |
88dd9ad4 | 351 | radiusMin = FromEmissionToCathode(); |
5cb4dfc3 | 352 | } |
353 | ||
354 | // Second value... | |
355 | ||
88dd9ad4 | 356 | FindThetaAtQuartz(thetaCerMax); |
5cb4dfc3 | 357 | if(GetThetaAtQuartz() == 999.) |
358 | { | |
88dd9ad4 | 359 | radiusMax = 999.; |
5cb4dfc3 | 360 | } |
361 | else | |
362 | { | |
363 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
88dd9ad4 | 364 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 365 | |
88dd9ad4 | 366 | radiusMax = FromEmissionToCathode(); |
5cb4dfc3 | 367 | } |
368 | // Mean value... | |
369 | ||
88dd9ad4 | 370 | thetaCerMean = (thetaCerMax + thetaCerMin)/2; |
5cb4dfc3 | 371 | |
88dd9ad4 | 372 | FindThetaAtQuartz(thetaCerMean); |
5cb4dfc3 | 373 | if(GetThetaAtQuartz() == 999.) |
374 | { | |
88dd9ad4 | 375 | radiusMean = 999.; |
5cb4dfc3 | 376 | } |
377 | else | |
378 | { | |
379 | SetThetaPhotonInDRS(GetThetaAtQuartz()); | |
88dd9ad4 | 380 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 381 | |
88dd9ad4 | 382 | radiusMean = FromEmissionToCathode(); |
5cb4dfc3 | 383 | } |
384 | ||
101624cd | 385 | // AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax)); |
5cb4dfc3 | 386 | |
88dd9ad4 | 387 | while (TMath::Abs(radiusMean-distPoint) > kTollerance) |
5cb4dfc3 | 388 | { |
389 | ||
88dd9ad4 | 390 | if((radiusMin-distPoint)*(radiusMean-distPoint) < 0) thetaCerMax = thetaCerMean; |
391 | if((radiusMin-distPoint)*(radiusMean-distPoint) > 0) { | |
5cb4dfc3 | 392 | |
88dd9ad4 | 393 | thetaCerMin = thetaCerMean; |
5cb4dfc3 | 394 | |
88dd9ad4 | 395 | FindThetaAtQuartz(thetaCerMin); |
5cb4dfc3 | 396 | SetThetaPhotonInDRS(GetThetaAtQuartz()); |
88dd9ad4 | 397 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 398 | |
88dd9ad4 | 399 | radiusMin =FromEmissionToCathode(); |
5cb4dfc3 | 400 | } |
401 | ||
88dd9ad4 | 402 | thetaCerMean = (thetaCerMax + thetaCerMin)/2; |
5cb4dfc3 | 403 | |
88dd9ad4 | 404 | FindThetaAtQuartz(thetaCerMean); |
5cb4dfc3 | 405 | SetThetaPhotonInDRS(GetThetaAtQuartz()); |
88dd9ad4 | 406 | SetPhiPhotonInDRS(phiPoint); |
5cb4dfc3 | 407 | |
88dd9ad4 | 408 | radiusMean = FromEmissionToCathode(); |
5cb4dfc3 | 409 | |
410 | nIteration++; | |
411 | if(nIteration>=50) { | |
d0831219 | 412 | // AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration)); |
5cb4dfc3 | 413 | SetThetaPhotonCerenkov(999.); |
414 | return; | |
415 | } | |
416 | } | |
417 | ||
101624cd | 418 | // AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean)); |
88dd9ad4 | 419 | SetThetaPhotonCerenkov(thetaCerMean); |
5cb4dfc3 | 420 | |
421 | } | |
101624cd | 422 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 423 | void AliRICHRecon::FindAreaAndPortionOfRing() |
424 | { | |
88dd9ad4 | 425 | //find fraction of the ring accepted by the RICH |
5cb4dfc3 | 426 | |
88dd9ad4 | 427 | Float_t xPoint[NPointsOfRing], yPoint[NPointsOfRing]; |
5cb4dfc3 | 428 | |
88dd9ad4 | 429 | // Float_t xtoentr = GetEntranceX(); |
430 | // Float_t ytoentr = GetEntranceY(); | |
431 | Float_t shiftX = GetShiftX(); | |
432 | Float_t shiftY = GetShiftY(); | |
5cb4dfc3 | 433 | |
88dd9ad4 | 434 | Float_t xemiss = GetXCoordOfEmission(); |
435 | Float_t yemiss = GetYCoordOfEmission(); | |
5cb4dfc3 | 436 | |
88dd9ad4 | 437 | Float_t x0 = xemiss + shiftX; |
438 | Float_t y0 = yemiss + shiftY; | |
5cb4dfc3 | 439 | |
5cb4dfc3 | 440 | |
5cb4dfc3 | 441 | |
db910db9 | 442 | SetEmissionPoint(AliRICHParam::RadThick()/2.); |
5cb4dfc3 | 443 | |
88dd9ad4 | 444 | Float_t theta = GetThetaOfRing(); |
5cb4dfc3 | 445 | |
446 | Int_t nPoints = 0; | |
88dd9ad4 | 447 | Int_t nPsiAccepted = 0; |
448 | Int_t nPsiTotal = 0; | |
5cb4dfc3 | 449 | |
db910db9 | 450 | for(Int_t i=0;i<NPointsOfRing-1;i++){ |
88dd9ad4 | 451 | Float_t psi = 2*TMath::Pi()*i/NPointsOfRing; |
5cb4dfc3 | 452 | |
88dd9ad4 | 453 | SetThetaPhotonInTRS(theta); |
454 | SetPhiPhotonInTRS(psi); | |
5cb4dfc3 | 455 | FindPhotonAnglesInDRS(); |
456 | ||
88dd9ad4 | 457 | Float_t radius = FromEmissionToCathode(); |
458 | if (radius == 999.) continue; | |
5cb4dfc3 | 459 | |
88dd9ad4 | 460 | nPsiTotal++; |
5cb4dfc3 | 461 | |
88dd9ad4 | 462 | Float_t xPointRing = GetXPointOnCathode() + shiftX; |
463 | Float_t yPointRing = GetYPointOnCathode() + shiftY; | |
5cb4dfc3 | 464 | |
88dd9ad4 | 465 | SetDetectorWhereX(xPointRing); |
466 | SetDetectorWhereY(yPointRing); | |
5cb4dfc3 | 467 | |
88dd9ad4 | 468 | Int_t zone = CheckDetectorAcceptance(); |
5cb4dfc3 | 469 | |
910735ae | 470 | AliDebug(1,Form("acceptance to detector zone -> %d",zone)); |
5cb4dfc3 | 471 | |
db910db9 | 472 | if (zone != 0){ |
473 | FindIntersectionWithDetector(); | |
474 | xPoint[nPoints] = GetIntersectionX(); yPoint[nPoints] = GetIntersectionY(); | |
475 | }else{ | |
476 | xPoint[nPoints] = xPointRing; yPoint[nPoints] = yPointRing; | |
477 | nPsiAccepted++; | |
478 | } | |
5cb4dfc3 | 479 | nPoints++; |
db910db9 | 480 | } |
5cb4dfc3 | 481 | |
db910db9 | 482 | xPoint[nPoints] = xPoint[0]; yPoint[nPoints] = yPoint[0]; |
5cb4dfc3 | 483 | |
484 | // find area... | |
485 | ||
88dd9ad4 | 486 | Float_t area = 0; |
5cb4dfc3 | 487 | |
488 | for (Int_t i = 0; i < nPoints; i++) | |
489 | { | |
88dd9ad4 | 490 | area += TMath::Abs((xPoint[i]-x0)*(yPoint[i+1]-y0) - (xPoint[i+1]-x0)*(yPoint[i]-y0)); |
5cb4dfc3 | 491 | } |
492 | ||
88dd9ad4 | 493 | area *= 0.5; |
5cb4dfc3 | 494 | |
adcdf4dd | 495 | Float_t portionOfRing = 0; |
496 | if (nPsiTotal>0) | |
497 | portionOfRing = ((Float_t)nPsiAccepted)/((Float_t)(nPsiTotal)); | |
5cb4dfc3 | 498 | |
5cb4dfc3 | 499 | |
88dd9ad4 | 500 | SetAreaOfRing(area); |
501 | SetPortionOfRing(portionOfRing); | |
db910db9 | 502 | }//FindAreaAndPortionOfRing() |
101624cd | 503 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 504 | void AliRICHRecon::FindIntersectionWithDetector() |
505 | { | |
88dd9ad4 | 506 | // find ring intersection with CsI edges |
5cb4dfc3 | 507 | |
88dd9ad4 | 508 | Float_t xIntersect, yIntersect; |
5cb4dfc3 | 509 | Float_t x1, x2, y1, y2; |
510 | ||
88dd9ad4 | 511 | Float_t shiftX = GetShiftX(); |
512 | Float_t shiftY = GetShiftY(); | |
5cb4dfc3 | 513 | |
88dd9ad4 | 514 | Float_t xPoint = GetXPointOnCathode() + shiftX; |
515 | Float_t yPoint = GetYPointOnCathode() + shiftY; | |
5cb4dfc3 | 516 | |
88dd9ad4 | 517 | Float_t xemiss = GetXCoordOfEmission(); |
518 | Float_t yemiss = GetYCoordOfEmission(); | |
5cb4dfc3 | 519 | |
88dd9ad4 | 520 | Float_t phi = GetPhiPhotonInDRS(); |
521 | Float_t m = tan(phi); | |
5cb4dfc3 | 522 | |
88dd9ad4 | 523 | Float_t x0 = xemiss + shiftX; |
524 | Float_t y0 = yemiss + shiftY; | |
5cb4dfc3 | 525 | |
88dd9ad4 | 526 | if(xPoint > x0) |
5cb4dfc3 | 527 | { |
528 | x1 = x0; | |
88dd9ad4 | 529 | x2 = xPoint; |
5cb4dfc3 | 530 | } |
531 | else | |
532 | { | |
533 | x2 = x0; | |
88dd9ad4 | 534 | x1 = xPoint; |
5cb4dfc3 | 535 | } |
88dd9ad4 | 536 | if(yPoint > y0) |
5cb4dfc3 | 537 | { |
538 | y1 = y0; | |
88dd9ad4 | 539 | y2 = yPoint; |
5cb4dfc3 | 540 | } |
541 | else | |
542 | { | |
543 | y2 = y0; | |
88dd9ad4 | 544 | y1 = yPoint; |
5cb4dfc3 | 545 | } |
546 | // | |
db910db9 | 547 | xIntersect = AliRICHParam::PcSizeX(); |
88dd9ad4 | 548 | yIntersect = m*(xIntersect - x0) + y0; |
db910db9 | 549 | if (yIntersect >= 0 && yIntersect <= AliRICHParam::PcSizeY() && xIntersect >= x1 && xIntersect <= x2) |
5cb4dfc3 | 550 | { |
88dd9ad4 | 551 | SetIntersectionX(xIntersect); |
552 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 553 | return; |
554 | } | |
555 | // | |
db910db9 | 556 | xIntersect = 0; |
88dd9ad4 | 557 | yIntersect = m*(xIntersect - x0) + y0; |
db910db9 | 558 | if (yIntersect >= 0 && yIntersect <= AliRICHParam::PcSizeY() && xIntersect >= x1 && xIntersect <= x2) |
5cb4dfc3 | 559 | { |
88dd9ad4 | 560 | SetIntersectionX(xIntersect); |
561 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 562 | return; |
563 | } | |
564 | // | |
db910db9 | 565 | yIntersect = AliRICHParam::PcSizeY(); |
88dd9ad4 | 566 | xIntersect = (yIntersect - y0)/m + x0; |
db910db9 | 567 | if (xIntersect >= 0 && xIntersect <= AliRICHParam::PcSizeX() && yIntersect >= y1 && yIntersect <= y2) |
5cb4dfc3 | 568 | { |
88dd9ad4 | 569 | SetIntersectionX(xIntersect); |
570 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 571 | return; |
572 | } | |
573 | // | |
db910db9 | 574 | yIntersect = 0; |
88dd9ad4 | 575 | xIntersect = (yIntersect - y0)/m + x0; |
db910db9 | 576 | if (xIntersect >= 0 && xIntersect <= AliRICHParam::PcSizeX() && yIntersect >= y1 && yIntersect <= y2) |
5cb4dfc3 | 577 | { |
88dd9ad4 | 578 | SetIntersectionX(xIntersect); |
579 | SetIntersectionY(yIntersect); | |
5cb4dfc3 | 580 | return; |
581 | } | |
5cb4dfc3 | 582 | } |
db910db9 | 583 | |
b068561d | 584 | //__________________________________________________________________________________________________ |
88dd9ad4 | 585 | Int_t AliRICHRecon::CheckDetectorAcceptance() const |
5cb4dfc3 | 586 | { |
88dd9ad4 | 587 | // check for the acceptance |
5cb4dfc3 | 588 | |
589 | // crosses X -2.6 2.6 cm | |
590 | // crosses Y -1 1 cm | |
591 | ||
88dd9ad4 | 592 | Float_t xcoord = GetDetectorWhereX(); |
593 | Float_t ycoord = GetDetectorWhereY(); | |
5cb4dfc3 | 594 | |
db910db9 | 595 | if(xcoord > AliRICHParam::PcSizeX()) |
5cb4dfc3 | 596 | { |
db910db9 | 597 | if(ycoord > AliRICHParam::PcSizeY()) return 2; |
598 | if(ycoord > 0 && ycoord < AliRICHParam::PcSizeY()) return 3; | |
599 | if(ycoord < 0) return 4; | |
5cb4dfc3 | 600 | } |
db910db9 | 601 | if(xcoord < 0) |
5cb4dfc3 | 602 | { |
db910db9 | 603 | if(ycoord > AliRICHParam::PcSizeY()) return 8; |
604 | if(ycoord > 0 && ycoord < AliRICHParam::PcSizeY()) return 7; | |
605 | if(ycoord < 0) return 6; | |
5cb4dfc3 | 606 | } |
db910db9 | 607 | if(xcoord > 0 && xcoord < AliRICHParam::PcSizeX()) |
5cb4dfc3 | 608 | { |
db910db9 | 609 | if(ycoord > AliRICHParam::PcSizeY()) return 1; |
610 | if(ycoord > 0 && ycoord < AliRICHParam::PcSizeY()) return 0; | |
611 | if(ycoord < 0) return 5; | |
5cb4dfc3 | 612 | } |
613 | return 999; | |
614 | } | |
b068561d | 615 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 616 | void AliRICHRecon::FindPhotonAnglesInDRS() |
617 | { | |
618 | // Setup the rotation matrix of the track... | |
619 | ||
88dd9ad4 | 620 | TRotation mtheta; |
621 | TRotation mphi; | |
622 | TRotation minv; | |
623 | TRotation mrot; | |
5cb4dfc3 | 624 | |
88dd9ad4 | 625 | Float_t trackTheta = GetTrackTheta(); |
626 | Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 627 | |
88dd9ad4 | 628 | mtheta.RotateY(trackTheta); |
629 | mphi.RotateZ(trackPhi); | |
5cb4dfc3 | 630 | |
88dd9ad4 | 631 | mrot = mphi * mtheta; |
632 | // minv = mrot.Inverse(); | |
5cb4dfc3 | 633 | |
88dd9ad4 | 634 | TVector3 photonInRadiator(1,1,1); |
5cb4dfc3 | 635 | |
88dd9ad4 | 636 | Float_t thetaCerenkov = GetThetaPhotonInTRS(); |
637 | Float_t phiCerenkov = GetPhiPhotonInTRS(); | |
5cb4dfc3 | 638 | |
88dd9ad4 | 639 | photonInRadiator.SetTheta(thetaCerenkov); |
640 | photonInRadiator.SetPhi(phiCerenkov); | |
641 | photonInRadiator = mrot * photonInRadiator; | |
642 | Float_t theta = photonInRadiator.Theta(); | |
643 | Float_t phi = photonInRadiator.Phi(); | |
644 | SetThetaPhotonInDRS(theta); | |
645 | SetPhiPhotonInDRS(phi); | |
5cb4dfc3 | 646 | |
647 | } | |
101624cd | 648 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 649 | Float_t AliRICHRecon::FromEmissionToCathode() |
650 | { | |
db910db9 | 651 | // Trace current photon from emission point somewhere in radiator to photocathode |
652 | // Arguments: none | |
653 | // Returns: | |
5cb4dfc3 | 654 | |
655 | Float_t nfreon, nquartz, ngas; | |
656 | ||
4c3a1f92 | 657 | nfreon = AliRICHParam::Instance()->IdxC6F14(AliRICHParam::Instance()->EckovMean()); |
658 | nquartz = AliRICHParam::Instance()->IdxSiO2(AliRICHParam::Instance()->EckovMean()); | |
659 | ngas = AliRICHParam::Instance()->IdxCH4(AliRICHParam::Instance()->EckovMean()); | |
5cb4dfc3 | 660 | |
88dd9ad4 | 661 | Float_t trackTheta = GetTrackTheta(); |
662 | Float_t trackPhi = GetTrackPhi(); | |
663 | Float_t lengthOfEmissionPoint = GetEmissionPoint(); | |
5cb4dfc3 | 664 | |
88dd9ad4 | 665 | Float_t theta = GetThetaPhotonInDRS(); |
666 | Float_t phi = GetPhiPhotonInDRS(); | |
5cb4dfc3 | 667 | |
88dd9ad4 | 668 | Float_t xemiss = lengthOfEmissionPoint*tan(trackTheta)*cos(trackPhi); |
669 | Float_t yemiss = lengthOfEmissionPoint*tan(trackTheta)*sin(trackPhi); | |
5cb4dfc3 | 670 | |
88dd9ad4 | 671 | SetXCoordOfEmission(xemiss); |
672 | SetYCoordOfEmission(yemiss); | |
5cb4dfc3 | 673 | |
88dd9ad4 | 674 | Float_t thetaquar = SnellAngle(nfreon, nquartz, theta); |
5cb4dfc3 | 675 | |
676 | if(thetaquar == 999.) | |
677 | { | |
678 | SetXPointOnCathode(999.); | |
679 | SetYPointOnCathode(999.); | |
680 | return thetaquar; | |
681 | } | |
682 | ||
683 | Float_t thetagap = SnellAngle( nquartz, ngas, thetaquar); | |
684 | ||
685 | if(thetagap == 999.) | |
686 | { | |
687 | SetXPointOnCathode(999.); | |
688 | SetYPointOnCathode(999.); | |
689 | return thetagap; | |
690 | } | |
691 | ||
db910db9 | 692 | Float_t xw = (AliRICHParam::RadThick() - lengthOfEmissionPoint)*cos(phi)*tan(theta); |
693 | Float_t xq = AliRICHParam::WinThick()*cos(phi)*tan(thetaquar); | |
694 | Float_t xg = AliRICHParam::Pc2Win()*cos(phi)*tan(thetagap); | |
695 | Float_t yw = (AliRICHParam::RadThick() - lengthOfEmissionPoint)*sin(phi)*tan(theta); | |
696 | Float_t yq = AliRICHParam::WinThick()*sin(phi)*tan(thetaquar); | |
697 | Float_t yg = AliRICHParam::Pc2Win()*sin(phi)*tan(thetagap); | |
5cb4dfc3 | 698 | |
5cb4dfc3 | 699 | |
88dd9ad4 | 700 | Float_t xtot = xemiss + xw + xq + xg; |
701 | Float_t ytot = yemiss + yw + yq + yg; | |
5cb4dfc3 | 702 | |
703 | SetXPointOnCathode(xtot); | |
704 | SetYPointOnCathode(ytot); | |
705 | ||
5cb4dfc3 | 706 | |
db910db9 | 707 | Float_t distanceFromEntrance = TMath::Sqrt(TMath::Power(fPhotonLimitX,2)+TMath::Power(fPhotonLimitY,2)); |
5cb4dfc3 | 708 | |
88dd9ad4 | 709 | return distanceFromEntrance; |
5cb4dfc3 | 710 | |
711 | } | |
101624cd | 712 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 713 | void AliRICHRecon::FindPhiPoint() |
714 | { | |
88dd9ad4 | 715 | //find phi of generated point |
5cb4dfc3 | 716 | |
88dd9ad4 | 717 | Float_t xtoentr = GetEntranceX(); |
718 | Float_t ytoentr = GetEntranceY(); | |
5cb4dfc3 | 719 | |
88dd9ad4 | 720 | Float_t trackTheta = GetTrackTheta(); |
721 | Float_t trackPhi = GetTrackPhi(); | |
5cb4dfc3 | 722 | |
88dd9ad4 | 723 | Float_t emissionPoint = GetEmissionPoint(); |
5cb4dfc3 | 724 | |
88dd9ad4 | 725 | Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi); |
726 | Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi); | |
101624cd | 727 | Float_t phi = atan2(argY,argX); |
5cb4dfc3 | 728 | |
998b831f | 729 | SetPhiPoint(phi); |
5cb4dfc3 | 730 | |
731 | } | |
101624cd | 732 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 733 | Float_t AliRICHRecon::Cerenkovangle(Float_t n, Float_t beta) |
734 | { | |
88dd9ad4 | 735 | // cerenkov angle from n and beta |
5cb4dfc3 | 736 | |
737 | // Compute the cerenkov angle | |
738 | ||
739 | Float_t thetacer; | |
740 | ||
741 | if((n*beta)<1.) { | |
742 | thetacer = 999.; | |
743 | // cout << " warning in Cerenkoangle !!!!!! " << endl; | |
744 | return thetacer; | |
745 | } | |
746 | ||
747 | thetacer = acos (1./(n*beta)); | |
748 | return thetacer; | |
749 | } | |
101624cd | 750 | //__________________________________________________________________________________________________ |
5cb4dfc3 | 751 | Float_t AliRICHRecon::SnellAngle(Float_t n1, Float_t n2, Float_t theta1) |
88dd9ad4 | 752 | { |
753 | // Snell law | |
5cb4dfc3 | 754 | |
755 | // Compute the Snell angle | |
756 | ||
757 | Float_t sinrefractangle; | |
758 | Float_t refractangle; | |
759 | ||
760 | sinrefractangle = (n1/n2)*sin(theta1); | |
761 | ||
762 | if(sinrefractangle>1.) { | |
763 | // cout << " PROBLEMS IN SNELL ANGLE !!!!! " << endl; | |
764 | refractangle = 999.; | |
765 | return refractangle; | |
766 | } | |
767 | ||
768 | refractangle = asin(sinrefractangle); | |
769 | return refractangle; | |
770 | } | |
101624cd | 771 | //__________________________________________________________________________________________________ |
db910db9 | 772 | Double_t AliRICHRecon::HoughResponse() |
88dd9ad4 | 773 | { |
db910db9 | 774 | // |
775 | // | |
776 | // | |
fab9e039 | 777 | Int_t nChannels = (Int_t)(fThetaMax/fDTheta+0.5); |
db910db9 | 778 | TH1F *phots = new TH1F("phots" ,"phots" ,nChannels,0,fThetaMax); |
779 | TH1F *photsw = new TH1F("photsw" ,"photsw" ,nChannels,0,fThetaMax); | |
780 | TH1F *resultw = new TH1F("resultw","resultw",nChannels,0,fThetaMax); | |
910735ae | 781 | Int_t nBin = (Int_t)(fThetaMax/fDTheta); |
782 | Int_t nCorrBand = (Int_t)(fWindowWidth/(2*fDTheta)); | |
783 | AliDebug(1,Form("Ring reconstruction for track with theta %f",GetTrackTheta()*TMath::RadToDeg())); | |
db910db9 | 784 | for (Int_t kPhot=0; kPhot< GetPhotonsNumber(); kPhot++){ |
910735ae | 785 | SetPhotonIndex(kPhot); |
786 | Double_t angle = GetPhotonEta(); | |
fab9e039 | 787 | if(angle<0||angle>fThetaMax) continue; |
910735ae | 788 | phots->Fill(angle); |
789 | Int_t bin = (Int_t)(0.5+angle/(fDTheta)); | |
790 | Double_t weight=1.; | |
791 | if(fIsWEIGHT){ | |
792 | Double_t lowerlimit = ((Float_t)bin)*fDTheta - 0.5*fDTheta; | |
793 | SetThetaOfRing(lowerlimit); | |
794 | FindAreaAndPortionOfRing(); | |
795 | Float_t area1 = GetAreaOfRing(); | |
796 | Double_t upperlimit = ((Float_t)bin)*fDTheta + 0.5*fDTheta; | |
797 | SetThetaOfRing(upperlimit); | |
798 | FindAreaAndPortionOfRing(); | |
799 | Float_t area2 = GetAreaOfRing(); | |
800 | AliDebug(1,Form("lowerlimit %f area %f ; upperlimit %f area %f",lowerlimit,area1,upperlimit,area2)); | |
801 | Float_t diffarea = area2 - area1; | |
802 | if(diffarea>0){weight = 1./(area2-area1);}else{weight = 1.;} | |
803 | } | |
fab9e039 | 804 | AliDebug(1,Form("Calculated weight %f",weight)); |
910735ae | 805 | photsw->Fill(angle,weight); |
806 | SetPhotonWeight(weight); | |
807 | } | |
808 | for (Int_t i=1; i<=nBin;i++){ | |
809 | Int_t bin1= i-nCorrBand; | |
810 | Int_t bin2= i+nCorrBand; | |
811 | if(bin1<1) bin1=1; | |
812 | if(bin2>nBin)bin2=nBin; | |
813 | Double_t sumPhots=phots->Integral(bin1,bin2); | |
814 | if(sumPhots<fMinNumPhots) continue; // cut on minimum n. of photons per ring | |
815 | Double_t sumPhotsw=photsw->Integral(bin1,bin2); | |
fab9e039 | 816 | resultw->Fill((Float_t)((i+0.5)*fDTheta),sumPhotsw); |
db910db9 | 817 | } |
818 | // evaluate the "BEST" theta ckov as the maximum value of histogramm | |
910735ae | 819 | Float_t *pVec = resultw->GetArray(); |
820 | Int_t locMax = TMath::LocMax(nBin,pVec); | |
db910db9 | 821 | phots->Delete();photsw->Delete();resultw->Delete(); // Reset and delete objects |
822 | ||
823 | return (Double_t)(locMax*fDTheta+0.5*fDTheta); //final most probable track theta ckov | |
910735ae | 824 | }//HoughResponse |
101624cd | 825 | //__________________________________________________________________________________________________ |
4c3a1f92 | 826 | void AliRICHRecon::FindThetaCerenkov(Int_t iNclus) |
5cb4dfc3 | 827 | { |
84093f6f | 828 | // Loops on all Ckov candidates and estimates the best Theta Ckov for a ring formed by those candidates. Also estimates an error for that Theat Ckov |
829 | // collecting errors for all single Ckov candidates thetas. (Assuming they are independent) | |
4c3a1f92 | 830 | // Arguments: iNclus- total number of clusters in chamber for background estimation |
84093f6f | 831 | // Return: none |
5cb4dfc3 | 832 | |
833 | Float_t wei = 0.; | |
88dd9ad4 | 834 | Float_t weightThetaCerenkov = 0.; |
5cb4dfc3 | 835 | |
fab9e039 | 836 | Double_t etaMin=9999.,etaMax=0.; |
84093f6f | 837 | Double_t sigma2 = 0; //to collect error squared for this ring |
838 | ||
db910db9 | 839 | for(Int_t i=0;i<GetPhotonsNumber();i++){ |
840 | SetPhotonIndex(i); | |
841 | if(GetPhotonFlag() == 2){ | |
84093f6f | 842 | Float_t photonEta = GetPhotonEta(); |
db910db9 | 843 | if(photonEta<etaMin) etaMin=photonEta; |
844 | if(photonEta>etaMax) etaMax=photonEta; | |
84093f6f | 845 | Float_t photonWeight = GetPhotonWeight(); |
846 | weightThetaCerenkov += photonEta*photonWeight; | |
847 | wei += photonWeight; | |
848 | //here comes sigma of the reconstructed ring | |
849 | ||
850 | //Double_t phiref=(GetPhiPoint()-GetTrackPhi()); | |
851 | if(GetPhotonEta()<=0) continue;//?????????????????Flag photos = 2 may imply CkovEta = 0?????????????? | |
852 | //??????????? Look at SetPhoton Flag method | |
853 | Double_t phiref=GetTrackPhi(); | |
854 | ||
4c3a1f92 | 855 | Double_t beta = 1./(TMath::Cos(GetPhotonEta())*AliRICHParam::Instance()->IdxC6F14(AliRICHParam::EckovMean())); |
84093f6f | 856 | sigma2 += 1./AliRICHParam::SigmaSinglePhotonFormula(GetPhotonEta(),GetPhiPoint(),GetTrackTheta(),phiref,beta); |
857 | } | |
db910db9 | 858 | } |
84093f6f | 859 | |
860 | if(sigma2>0) SetRingSigma2(1./sigma2); | |
861 | else SetRingSigma2(1e10); | |
862 | ||
fab9e039 | 863 | if(wei != 0.) weightThetaCerenkov /= wei; else weightThetaCerenkov = 0.; |
88dd9ad4 | 864 | SetThetaCerenkov(weightThetaCerenkov); |
5cb4dfc3 | 865 | |
fab9e039 | 866 | // estimate of the n. of bkg photons |
867 | SetThetaOfRing(etaMin); FindAreaAndPortionOfRing(); Double_t internalArea = GetAreaOfRing(); | |
868 | SetThetaOfRing(etaMax); FindAreaAndPortionOfRing(); Double_t externalArea = GetAreaOfRing(); | |
869 | ||
870 | Double_t effArea = (AliRICHParam::PcSizeX()-AliRICHParam::DeadZone())*(AliRICHParam::PcSizeY()-2*AliRICHParam::DeadZone()); | |
4c3a1f92 | 871 | Double_t nPhotBKG = (externalArea-internalArea)/effArea*iNclus;//????? is the division right? |
fab9e039 | 872 | if(nPhotBKG<0) nPhotBKG=0; //just protection from funny angles... |
873 | SetPhotBKG(nPhotBKG); | |
fab9e039 | 874 | |
101624cd | 875 | AliDebug(1,Form(" thetac weighted -> %f",weightThetaCerenkov)); |
4c3a1f92 | 876 | }//FindThetaCerenkov() |
db910db9 | 877 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
878 | Int_t AliRICHRecon::FlagPhotons(Double_t thetaCkovHough) | |
5cb4dfc3 | 879 | { |
db910db9 | 880 | // flag photon candidates if their individual theta ckov inside the window around theta ckov of Hough transform |
881 | // Arguments: thetaCkovHough- value of most probable theta ckov for track as returned by HoughResponse() | |
882 | // Returns: number of photon candidates happened to be inside the window | |
5cb4dfc3 | 883 | |
db910db9 | 884 | Int_t steps = (Int_t)((thetaCkovHough - fThetaMin)/ fDTheta); //how many times we need to have fDTheta to fill the distance betwee fThetaMin and thetaCkovHough |
5cb4dfc3 | 885 | |
886 | Float_t tmin = fThetaMin + (Float_t)(steps - 1)*fDTheta; | |
887 | Float_t tmax = fThetaMin + (Float_t)(steps)*fDTheta; | |
888 | Float_t tavg = 0.5*(tmin+tmax); | |
889 | ||
db910db9 | 890 | tmin = tavg - 0.5*fWindowWidth; tmax = tavg + 0.5*fWindowWidth; |
5cb4dfc3 | 891 | |
db910db9 | 892 | Int_t iInsideCnt = 0; //count photons which theta inside prdefined window |
893 | for(Int_t i=0;i<GetPhotonsNumber();i++){//photon candidates loop | |
894 | SetPhotonIndex(i); Float_t photonEta = GetPhotonEta(); | |
895 | if(photonEta == -999.) continue; | |
84093f6f | 896 | if(photonEta >= tmin && photonEta <= tmax) { |
897 | SetPhotonFlag(2); | |
898 | iInsideCnt++; | |
899 | } | |
db910db9 | 900 | } |
901 | return iInsideCnt; | |
902 | }//FlagPhotons |