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