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