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