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