Fixed forward declarations.
[u/mrichter/AliRoot.git] / RICH / AliRICHPatRec.cxx
CommitLineData
e7257cad 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 $Log$
488e98ba 18 Revision 1.4 2000/06/30 16:33:43 dibari
19 Several changes (ring drawing, fiducial selection, etc.)
20
e0b63a71 21 Revision 1.3 2000/06/15 15:47:12 jbarbosa
22 Corrected compilation errors on HP-UX (replaced pow with TMath::Power)
23
00df6e79 24 Revision 1.2 2000/06/12 15:26:09 jbarbosa
25 Cleaned up version.
26
237c933d 27 Revision 1.1 2000/06/09 14:53:01 jbarbosa
28 Bari's pattern recognition algorithm
29
e7257cad 30*/
31
237c933d 32#include "AliRICHHit.h"
33#include "AliRICHCerenkov.h"
34#include "AliRICHPadHit.h"
35#include "AliRICHDigit.h"
36#include "AliRICHRawCluster.h"
37#include "AliRICHRecHit.h"
e7257cad 38#include "AliRun.h"
39#include "AliDetector.h"
40#include "AliRICH.h"
41#include "AliRICHPoints.h"
237c933d 42#include "AliRICHSegmentation.h"
e7257cad 43#include "AliRICHPatRec.h"
44#include "AliRICH.h"
45#include "AliRICHConst.h"
46#include "AliRICHPoints.h"
47#include "AliConst.h"
488e98ba 48#include "AliRICHHitMap.h"
237c933d 49
50#include <TParticle.h>
51#include <TMath.h>
52#include <TRandom.h>
53#include <TCanvas.h>
54#include <TH2.h>
e7257cad 55
56
57ClassImp(AliRICHPatRec)
58//___________________________________________
59AliRICHPatRec::AliRICHPatRec() : TObject()
60{
237c933d 61 // Default constructor
62
e7257cad 63 //fChambers = 0;
64}
65//___________________________________________
66AliRICHPatRec::AliRICHPatRec(const char *name, const char *title)
67 : TObject()
68{
237c933d 69 //Constructor for Bari's pattern recogniton method object
e7257cad 70}
71
72void AliRICHPatRec::PatRec()
73{
74
237c933d 75// Pattern recognition algorithm
76
e7257cad 77 AliRICHChamber* iChamber;
78 AliRICHSegmentation* segmentation;
79
237c933d 80 Int_t ntracks, ndigits[kNCH];
e7257cad 81 Int_t itr, ich, i;
237c933d 82 Int_t goodPhotons;
e7257cad 83 Int_t x,y,q;
84 Float_t rx,ry;
85 Int_t nent,status;
e0b63a71 86 Int_t padsUsedX[100];
87 Int_t padsUsedY[100];
e7257cad 88
e0b63a71 89 Float_t rechit[7];
e7257cad 90
91 printf("PatRec started\n");
92
237c933d 93 AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH");
94 TTree *treeH = gAlice->TreeH();
e7257cad 95
237c933d 96 ntracks =(Int_t) treeH->GetEntries();
e7257cad 97 // ntracks = 1;
98 for (itr=0; itr<ntracks; itr++) {
99
100 status = TrackParam(itr,ich);
101 if(status==1) continue;
e0b63a71 102 //printf(" theta %f phi %f track \n",fTrackTheta,fTrackPhi);
e7257cad 103 // ring->Fill(fTrackLoc[0],fTrackLoc[1],100.);
104
237c933d 105 iChamber = &(pRICH->Chamber(ich));
e7257cad 106 segmentation=iChamber->GetSegmentationModel();
107
108 nent=(Int_t)gAlice->TreeD()->GetEntries();
109 gAlice->TreeD()->GetEvent(nent-1);
237c933d 110 TClonesArray *pDigitss = pRICH->DigitsAddress(ich);
111 ndigits[ich] = pDigitss->GetEntriesFast();
e0b63a71 112 printf("Digits in chamber %d: %d\n",ich,ndigits[ich]);
e7257cad 113 AliRICHDigit *padI = 0;
114
237c933d 115 goodPhotons = 0;
e7257cad 116
117 for (Int_t dig=0;dig<ndigits[ich];dig++) {
237c933d 118 padI=(AliRICHDigit*) pDigitss->UncheckedAt(dig);
e7257cad 119 x=padI->fPadX;
120 y=padI->fPadY;
121 q=padI->fSignal;
122 segmentation->GetPadCxy(x,y,rx,ry);
123
e0b63a71 124 //printf("Pad coordinates x:%d, Real coordinates x:%f\n",x,rx);
125 //printf("Pad coordinates y:%d, Real coordinates y:%f\n",y,ry);
126
e7257cad 127 fXpad = rx-fXshift;
128 fYpad = ry-fYshift;
129 fQpad = q;
130
e7257cad 131 fCerenkovAnglePad = PhotonCerenkovAngle();
132 if(fCerenkovAnglePad==-999) continue;
133
134 if(!PhotonInBand()) continue;
135
e0b63a71 136 Int_t xpad;
137 Int_t ypad;
138
139 segmentation->GetPadIxy(fXpad,fYpad,xpad,ypad);
e7257cad 140
e0b63a71 141 padsUsedX[goodPhotons]=xpad;
142 padsUsedY[goodPhotons]=ypad;
143
237c933d 144 goodPhotons++;
e0b63a71 145 fEtaPhotons[goodPhotons-1] = fCerenkovAnglePad;
e7257cad 146 }
237c933d 147 fNumEtaPhotons = goodPhotons;
e7257cad 148
149 BackgroundEstimation();
150
e7257cad 151 HoughResponse();
e0b63a71 152 //CerenkovRingDrawing();
e7257cad 153
154 rechit[0] = 0;
155 rechit[1] = 0;
156 rechit[2] = fThetaCerenkov;
e0b63a71 157 rechit[3] = fXshift + fTrackLoc[0];
158 rechit[4] = fYshift + fTrackLoc[1];
159 rechit[5] = fEmissPoint;
160 rechit[6] = goodPhotons;
e7257cad 161
e0b63a71 162 //printf("Center coordinates:%f %f\n",rechit[3],rechit[4]);
e7257cad 163
e0b63a71 164 pRICH->AddRecHit(ich,rechit,fEtaPhotons,padsUsedX,padsUsedY);
e7257cad 165
e7257cad 166 }
167
168 gAlice->TreeR()->Fill();
169 TClonesArray *fRec;
237c933d 170 for (i=0;i<kNCH;i++) {
171 fRec=pRICH->RecHitsAddress(i);
e7257cad 172 int ndig=fRec->GetEntriesFast();
173 printf ("Chamber %d, rings %d\n",i,ndig);
174 }
237c933d 175 pRICH->ResetRecHits();
e7257cad 176
e7257cad 177}
178
179
180Int_t AliRICHPatRec::TrackParam(Int_t itr, Int_t &ich)
181{
182 // Get Local coordinates of track impact
183
184 AliRICHChamber* iChamber;
185 AliRICHSegmentation* segmentation;
186
187 Float_t trackglob[3];
188 Float_t trackloc[3];
189 Float_t thetatr;
190 Float_t phitr;
191 Float_t iloss;
192 Float_t part;
193 Float_t pX, pY, pZ;
194
195 printf("Calling TrackParam\n");
196
197 gAlice->ResetHits();
237c933d 198 TTree *treeH = gAlice->TreeH();
199 treeH->GetEvent(itr);
e7257cad 200
237c933d 201 AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH");
202 AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
e7257cad 203 if(mHit==0) return 1;
204 ich = mHit->fChamber-1;
205 trackglob[0] = mHit->fX;
206 trackglob[1] = mHit->fY;
207 trackglob[2] = mHit->fZ;
208 pX = mHit->fMomX;
209 pY = mHit->fMomY;
210 pZ = mHit->fMomZ;
00df6e79 211 fTrackMom = sqrt(TMath::Power(pX,2)+TMath::Power(pY,2)+TMath::Power(pZ,2));
e0b63a71 212 thetatr = (mHit->fTheta)*(Float_t)kDegrad;
e7257cad 213 phitr = mHit->fPhi*(Float_t)kDegrad;
214 iloss = mHit->fLoss;
215 part = mHit->fParticle;
216
237c933d 217 iChamber = &(pRICH->Chamber(ich));
e7257cad 218 iChamber->GlobaltoLocal(trackglob,trackloc);
219
220 segmentation=iChamber->GetSegmentationModel();
221
222 // retrieve geometrical params
223
224 AliRICHGeometry* fGeometry=iChamber->GetGeometryModel();
225
226 fRw = fGeometry->GetFreonThickness();
227 fQw = fGeometry->GetQuartzThickness();
e0b63a71 228 fTgap = fGeometry->GetGapThickness();
229 Float_t radiatorToPads= fGeometry->GetRadiatorToPads();
230 //+ fGeometry->GetProximityGapThickness();
e7257cad 231
e0b63a71 232 //printf("Distance to pads. From geometry:%f, From calculations:%f\n",radiatorToPads,fRw + fQw + fTgap);
233
234 //Float_t apar = (fRw + fQw + fTgap)*tan(thetatr);
235 Float_t apar = radiatorToPads*tan(thetatr);
e7257cad 236 fTrackLoc[0] = apar*cos(phitr);
237 fTrackLoc[1] = apar*sin(phitr);
e0b63a71 238 //fTrackLoc[2] = fRw + fQw + fTgap;
239 fTrackLoc[2] = radiatorToPads;
e7257cad 240 fTrackTheta = thetatr;
241 fTrackPhi = phitr;
242
243 fXshift = trackloc[0] - fTrackLoc[0];
244 fYshift = trackloc[2] - fTrackLoc[1];
245
246 return 0;
247}
248
249Float_t AliRICHPatRec::EstimationAtLimits(Float_t lim, Float_t radius,
250 Float_t phiphot)
251{
237c933d 252
253// Estimation of emission point
254
e7257cad 255 Float_t nquartz = 1.585;
256 Float_t ngas = 1.;
257 Float_t nfreon = 1.295;
258 Float_t value;
259
260 // printf("Calling EstimationLimits\n");
261
262 Float_t apar = (fRw -fEmissPoint + fQw + fTgap)*tan(fTrackTheta);
263 Float_t b1 = (fRw-fEmissPoint)*tan(lim);
00df6e79 264 Float_t b2 = fQw / sqrt(TMath::Power(nquartz,2)-TMath::Power(nfreon*sin(lim),2));
265 Float_t b3 = fTgap / sqrt(TMath::Power(ngas,2)-TMath::Power(nfreon*sin(lim),2));
e7257cad 266 Float_t bpar = b1 + nfreon*sin(lim)*(b2+b3);
00df6e79 267 value = TMath::Power(radius,2)
268 -TMath::Power((apar*cos(fTrackPhi)-bpar*cos(phiphot)),2)
269 -TMath::Power((apar*sin(fTrackPhi)-bpar*sin(phiphot)),2);
e7257cad 270 return value;
271}
272
273
274Float_t AliRICHPatRec::PhotonCerenkovAngle()
275{
276 // Cherenkov pad angle reconstruction
277
278 Float_t radius;
279 Float_t cherMin = 0;
280 Float_t cherMax = 0.8;
281 Float_t phiphot;
282 Float_t eps = 0.0001;
283 Int_t niterEmiss = 0;
284 Int_t niterEmissMax = 0;
237c933d 285 Float_t x1,x2,x3=0,p1,p2,p3;
e7257cad 286 Float_t argY,argX;
287 Int_t niterFun;
288
289 // printf("Calling PhotonCerenkovAngle\n");
290
00df6e79 291 radius = sqrt(TMath::Power(fTrackLoc[0]-fXpad,2)+TMath::Power(fTrackLoc[1]-fYpad,2));
e7257cad 292 fEmissPoint = fRw/2.; //Start value of EmissionPoint
293
294 while(niterEmiss<=niterEmissMax) {
295
296 niterFun = 0;
297 argY = fYpad - fEmissPoint*tan(fTrackTheta)*sin(fTrackPhi);
298 argX = fXpad - fEmissPoint*tan(fTrackTheta)*cos(fTrackPhi);
299 phiphot = atan2(argY,argX);
300 p1 = EstimationAtLimits(cherMin,radius,phiphot);
301 p2 = EstimationAtLimits(cherMax,radius,phiphot);
302 if(p1*p2>0)
303 {
304 // printf("PhotonCerenkovAngle failed\n");
305 return -999;
306 }
307
308 //start to find the Cherenkov pad angle
309 x1 = cherMin;
310 x2 = cherMax;
311 x3 = (x1+x2)/2.;
312 p3 = EstimationAtLimits(x3,radius,phiphot);
313 while(TMath::Abs(p3)>eps){
314 if(p1*p3<0) x2 = x3;
315 if(p1*p3>0) {
316 x1 = x3;
317 p1 = EstimationAtLimits(x1,radius,phiphot);
318 }
319 x3 = (x1+x2)/2.;
320 p3 = EstimationAtLimits(x3,radius,phiphot);
321 niterFun++;
322
323 if(niterFun>=1000) {
324 // printf(" max iterations in PhotonCerenkovAngle\n");
325 return x3;
326 }
327 }
328 // printf("niterFun %i \n",niterFun);
329 niterEmiss++;
330 if (niterEmiss != niterEmissMax+1) EmissionPoint();
331 }
332 /*
333 printf(" phiphot %f fXpad %f fYpad %f fEmiss %f \n",
334 phiphot,fXpad,fYpad,fEmissPoint);
335 */
336
337 return x3;
338
339}
340
341
342void AliRICHPatRec::EmissionPoint()
343{
237c933d 344
345// Find emission point
346
347 Float_t absorbtionLength=7.83*fRw; //absorption length in the freon (cm)
e7257cad 348 // 7.83 = -1/ln(T0) where
349 // T0->Trasmission freon at 180nm = 0.88 (Eph=6.85eV)
237c933d 350 Float_t photonLength, photonLengthMin, photonLengthMax;
e7257cad 351
237c933d 352 photonLength=exp(-fRw/(absorbtionLength*cos(fCerenkovAnglePad)));
353 photonLengthMin=fRw*photonLength/(1.-photonLength);
354 photonLengthMax=absorbtionLength*cos(fCerenkovAnglePad);
355 fEmissPoint = fRw + photonLengthMin - photonLengthMax;
e7257cad 356
357}
358
359void AliRICHPatRec::PhotonSelection(Int_t track, Int_t &nphot, Float_t &thetamean)
360{
237c933d 361
362// not implemented yet
363
e7257cad 364 printf("Calling PhotonSelection\n");
365}
366
367void AliRICHPatRec::BackgroundEstimation()
368{
237c933d 369
370// estimate background noise
371
372 Float_t stepEta = 0.001;
373 Float_t etaMinBkg = 0.72;
374 Float_t etaMaxBkg = 0.75;
375 Float_t etaMin = 0.;
376 Float_t etaMax = 0.75;
e7257cad 377 Float_t ngas = 1.;
378 Float_t nfreon = 1.295;
379
237c933d 380 Float_t etaStepMin,etaStepMax,etaStepAvg;
e7257cad 381 Int_t i,ip,nstep;
237c933d 382 Int_t numPhotBkg, numPhotonStep;
383 Float_t funBkg,areaBkg,normBkg;
384 Float_t densityBkg,storeBkg,numStore;
385 Float_t thetaSig;
e7257cad 386
237c933d 387 numPhotBkg = 0;
388 areaBkg = 0.;
e7257cad 389
237c933d 390 nstep = (int)((etaMaxBkg-etaMinBkg)/stepEta);
e7257cad 391
392 for (i=0;i<fNumEtaPhotons;i++) {
393
237c933d 394 if(fEtaPhotons[i]>etaMinBkg && fEtaPhotons[i]<etaMaxBkg) {
395 numPhotBkg++;
e7257cad 396 }
397 }
237c933d 398 if (numPhotBkg == 0) {
e7257cad 399 for (i=0;i<fNumEtaPhotons;i++) {
400 fWeightPhotons[i] = 1.;
401 }
402 return;
403 }
404
237c933d 405 // printf(" numPhotBkg %i ",numPhotBkg);
e7257cad 406
407 for (i=0;i<nstep;i++) {
237c933d 408 etaStepMin = etaMinBkg + (Float_t)(i)*stepEta;
409 etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;
410 etaStepAvg = 0.5*(etaStepMax + etaStepMin);
e7257cad 411 /*
00df6e79 412 funBkg = tan(etaStepAvg)*TMath::Power((1.+TMath::Power(tan(etaStepAvg),2)),
237c933d 413 5.52)-7.803 + 22.02*tan(etaStepAvg);
e7257cad 414 */
237c933d 415 thetaSig = asin(nfreon/ngas*sin(etaStepAvg));
00df6e79 416 funBkg = tan(thetaSig)*(1.+TMath::Power(tan(thetaSig),2))*nfreon
237c933d 417 /ngas*cos(etaStepAvg)/cos(thetaSig);
418 areaBkg += stepEta*funBkg;
e7257cad 419 }
420
237c933d 421 densityBkg = 0.95*(Float_t)(numPhotBkg)/areaBkg;
422 // printf(" densityBkg %f \n",densityBkg);
e7257cad 423
237c933d 424 nstep = (int)((etaMax-etaMin)/stepEta);
425 storeBkg = 0.;
426 numStore = 0;
e7257cad 427 for (i=0;i<nstep;i++) {
237c933d 428 etaStepMin = etaMinBkg + (Float_t)(i)*stepEta;
429 etaStepMax = etaMinBkg + (Float_t)(i+1)*stepEta;
430 etaStepAvg = 0.5*(etaStepMax + etaStepMin);
e7257cad 431 /*
00df6e79 432 funBkg = tan(etaStepAvg)*TMath::Power((1.+TMath::Power(tan(etaStepAvg),2)),
237c933d 433 5.52)-7.803 + 22.02*tan(etaStepAvg);
e7257cad 434 */
435
237c933d 436 thetaSig = asin(nfreon/ngas*sin(etaStepAvg));
00df6e79 437 funBkg = tan(thetaSig)*(1.+TMath::Power(tan(thetaSig),2))*nfreon
237c933d 438 /ngas*cos(etaStepAvg)/cos(thetaSig);
e7257cad 439
237c933d 440 areaBkg = stepEta*funBkg;
441 normBkg = densityBkg*areaBkg;
442 numPhotonStep = 0;
e7257cad 443 for (ip=0;ip<fNumEtaPhotons;ip++) {
237c933d 444 if(fEtaPhotons[ip]>etaStepMin && fEtaPhotons[ip]<etaStepMax) {
445 numPhotonStep++;
e7257cad 446 }
447 }
237c933d 448 if (numPhotonStep == 0) {
449 storeBkg += normBkg;
450 numStore++;
451 if (numStore>50) {
452 numStore = 0;
453 storeBkg = 0.;
e7257cad 454 }
455 }
237c933d 456 if (numPhotonStep == 0) continue;
e7257cad 457 for (ip=0;ip<fNumEtaPhotons;ip++) {
237c933d 458 if(fEtaPhotons[ip]>etaStepMin && fEtaPhotons[ip]<etaStepMax) {
459 normBkg +=storeBkg;
460 storeBkg = 0;
461 numStore = 0;
462 fWeightPhotons[ip] = 1. - normBkg/(Float_t)(numPhotonStep);
e7257cad 463 /*
237c933d 464 printf(" normBkg %f numPhotonStep %i fW %f \n",
465 normBkg, numPhotonStep, fWeightPhotons[ip]);
e7257cad 466 */
467 if(fWeightPhotons[ip]<0) fWeightPhotons[ip] = 0.;
468 }
469 }
470 }
471}
472
473
474void AliRICHPatRec::FlagPhotons(Int_t track, Float_t theta)
475{
237c933d 476
477// not implemented yet
478
e7257cad 479 printf("Calling FlagPhotons\n");
480}
481
482
483//////////////////////////////////////////
484
485
486
487
488
489Int_t AliRICHPatRec::PhotonInBand()
490{
491 //0=label for parameters giving internal band ellipse
492 //1=label for parameters giving external band ellipse
493
237c933d 494 Float_t imp[2], mass[2], energy[2], beta[2];
495 Float_t emissPointLength[2];
496 Float_t e1, e2, f1, f2;
e7257cad 497 Float_t nfreon[2], nquartz[2];
498 Int_t times;
e0b63a71 499 Float_t pointsOnCathode[3];
e7257cad 500
501 Float_t phpad, thetacer[2];
502 Float_t bandradius[2], padradius;
503
504 imp[0] = 5.0; //threshold momentum for the proton Cherenkov emission
505 imp[1] = 1.2;
506
507 mass[0] = 0.938; //proton mass
508 mass[1] = 0.139; //pion mass
509
237c933d 510 emissPointLength[0] = fRw-0.0001; //at the beginning of the radiator
511 emissPointLength[1] = 0.;//at the end of radiator
e7257cad 512
513 //parameters to calculate freon window refractive index vs. energy
514 Float_t a = 1.177;
515 Float_t b = 0.0172;
516
517 //parameters to calculate quartz window refractive index vs. energy
518 /*
519 Energ[0] = 5.6;
520 Energ[1] = 7.7;
521 */
237c933d 522 energy[0] = 5.0;
523 energy[1] = 8.0;
524 e1 = 10.666;
525 e2 = 18.125;
526 f1 = 46.411;
527 f2 = 228.71;
e7257cad 528
e7257cad 529 phpad = PhiPad();
530
531 for (times=0; times<=1; times++) {
532
237c933d 533 nfreon[times] = a+b*energy[times];
e0b63a71 534 //nfreon[times] = 1;
e7257cad 535
00df6e79 536 nquartz[times] = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(energy[times],2)))+
537 (f2/(TMath::Power(e2,2)-TMath::Power(energy[times],2))));
e7257cad 538
00df6e79 539 beta[times] = imp[times]/sqrt(TMath::Power(imp[times],2)+TMath::Power(mass[times],2));
e7257cad 540
541 thetacer[times] = CherenkovAngle( nfreon[times], beta[times]);
542
543 bandradius[times] = DistanceFromMip( nfreon[times], nquartz[times],
e0b63a71 544 emissPointLength[times],
545 thetacer[times], phpad, pointsOnCathode);
546 //printf(" ppp %f %f %f \n",pointsOnCathode);
547}
e7257cad 548
549 bandradius[0] -= 1.6;
550 bandradius[1] += 1.6;
00df6e79 551 padradius = sqrt(TMath::Power(fXpad,2)+TMath::Power(fYpad,2));
e7257cad 552 // printf(" rmin %f r %f rmax %f \n",bandradius[0],padradius,bandradius[1]);
553
554 if(padradius>=bandradius[0] && padradius<=bandradius[1]) return 1;
555 return 0;
556}
557
558Float_t AliRICHPatRec::DistanceFromMip(Float_t nfreon, Float_t nquartz,
237c933d 559 Float_t emissPointLength, Float_t thetacer,
e0b63a71 560 Float_t phpad, Float_t pointsOnCathode[3])
e7257cad 561{
e7257cad 562
237c933d 563// Find the distance to MIP impact
564
565 Float_t distanceValue;
566
567 TVector3 radExitPhot(1,1,1);//photon impact at the radiator exit with respect
e7257cad 568 //to local reference sistem with the origin in the MIP entrance
569
237c933d 570 TVector3 vectEmissPointLength(1,1,1);
571 Float_t magEmissPointLenght;
e7257cad 572
237c933d 573 TVector3 radExitPhot2(1,1,1);//photon impact at the radiator exit with respect
574 Float_t magRadExitPhot2;
e7257cad 575 //to a reference sistem with origin in the photon emission point and
576 //axes parallel to the MIP reference sistem
577
237c933d 578 TVector3 quarExitPhot(1,1,1);//photon impact at the quartz exit with respect
579 Float_t magQuarExitPhot;
e7257cad 580 //
237c933d 581 TVector3 gapExitPhot(1,1,1) ;
582 Float_t magGapExitPhot;
e7257cad 583 //
e0b63a71 584 TVector3 PhotocatExitPhot(1,1,1);
e7257cad 585 Double_t theta2;
586 Double_t thetarad , phirad ;
587 Double_t thetaquar, phiquar;
588 Double_t thetagap , phigap ;
589
590 Float_t ngas = 1.;
591
237c933d 592 magEmissPointLenght = emissPointLength/cos(fTrackTheta);
e7257cad 593
237c933d 594 vectEmissPointLength.SetMag(magEmissPointLenght);
595 vectEmissPointLength.SetTheta(fTrackTheta);
596 vectEmissPointLength.SetPhi(fTrackPhi);
e7257cad 597
598
237c933d 599 radExitPhot2.SetTheta(thetacer);
600 radExitPhot2.SetPhi(phpad);
e7257cad 601
602
603 TRotation r1;
604 TRotation r2;
605 TRotation r;
606
607 r1. RotateY(fTrackTheta);
608 r2. RotateZ(fTrackPhi);
609
610
611
612 r = r2 * r1;//rotation about the y axis by MIP theta incidence angle
613 //following by a rotation about the z axis by MIP phi incidence angle;
614
615
237c933d 616 radExitPhot2 = r * radExitPhot2;
617 theta2 = radExitPhot2.Theta();
618 magRadExitPhot2 = (fRw - vectEmissPointLength(2))/cos(theta2);
619 radExitPhot2.SetMag(magRadExitPhot2);
e7257cad 620
621
237c933d 622 radExitPhot = vectEmissPointLength + radExitPhot2;
623 thetarad = radExitPhot.Theta();
237c933d 624 phirad = radExitPhot.Phi(); //check on the original file //
e7257cad 625
626 thetaquar = SnellAngle( nfreon, nquartz, theta2);
237c933d 627 phiquar = radExitPhot2.Phi();
e7257cad 628 if(thetaquar == 999.) return thetaquar;
237c933d 629 magQuarExitPhot = fQw/cos(thetaquar);
630 quarExitPhot.SetMag( magQuarExitPhot);
631 quarExitPhot.SetTheta(thetaquar);
632 quarExitPhot.SetPhi(phiquar);
e7257cad 633
634 thetagap = SnellAngle( nquartz, ngas, thetaquar);
635 phigap = phiquar;
636 if(thetagap == 999.) return thetagap;
237c933d 637 magGapExitPhot = fTgap/cos(thetagap);
638 gapExitPhot.SetMag( magGapExitPhot);
639 gapExitPhot.SetTheta(thetagap);
640 gapExitPhot.SetPhi(phigap);
e7257cad 641
e0b63a71 642 PhotocatExitPhot = radExitPhot + quarExitPhot + gapExitPhot;
e7257cad 643
e0b63a71 644 distanceValue = sqrt(TMath::Power(PhotocatExitPhot(0),2)
645 +TMath::Power(PhotocatExitPhot(1),2));
646 pointsOnCathode[0] = (Float_t) PhotocatExitPhot(0) + fXshift - fTrackLoc[0];
647 pointsOnCathode[1] = (Float_t) PhotocatExitPhot(1) + fYshift - fTrackLoc[1];
648 pointsOnCathode[2] = (Float_t) PhotocatExitPhot(2);
649
650 //printf(" point in Distance.2. %f %f %f \n",pointsOnCathode[0],pointsOnCathode[1],pointsOnCathode[2]);
651
652 return distanceValue;
653
654 }
e7257cad 655
656Float_t AliRICHPatRec::PhiPad()
657{
237c933d 658
659// ??
660
e7257cad 661 Float_t zpad;
662 Float_t thetapad, phipad;
663 Float_t thetarot, phirot;
664
665 zpad = fRw + fQw + fTgap;
666
237c933d 667 TVector3 photonPad(fXpad, fYpad, zpad);
668 thetapad = photonPad.Theta();
669 phipad = photonPad.Phi();
e7257cad 670
671 TRotation r1;
672 TRotation r2;
673 TRotation r;
674
675 thetarot = - fTrackTheta;
676 phirot = - fTrackPhi;
677 r1. RotateZ(phirot);
678 r2. RotateY(thetarot);
679
680 r = r2 * r1;//rotation about the z axis by MIP -phi incidence angle
681 //following by a rotation about the y axis by MIP -theta incidence angle;
682
237c933d 683 photonPad = r * photonPad;
e7257cad 684
237c933d 685 phipad = photonPad.Phi();
e7257cad 686
687 return phipad;
688}
689
690Float_t AliRICHPatRec:: SnellAngle(Float_t n1, Float_t n2, Float_t theta1)
691{
237c933d 692
693// Compute the Snell angle
694
e7257cad 695 Float_t sinrefractangle;
696 Float_t refractangle;
697
698 sinrefractangle = (n1/n2)*sin(theta1);
699
700 if(sinrefractangle>1.) {
701 refractangle = 999.;
702 return refractangle;
703 }
704
705 refractangle = asin(sinrefractangle);
706 return refractangle;
707}
708
709Float_t AliRICHPatRec::CherenkovAngle(Float_t n, Float_t beta)
710{
237c933d 711
712// Compute the cerenkov angle
713
e7257cad 714 Float_t thetacer;
715
716 if((n*beta)<1.) {
717 thetacer = 999.;
718 return thetacer;
719 }
720
721 thetacer = acos (1./(n*beta));
722 return thetacer;
723}
724
725Float_t AliRICHPatRec::BetaCerenkov(Float_t n, Float_t theta)
726{
237c933d 727
728// Find beta
729
e7257cad 730 Float_t beta;
731
732 beta = 1./(n*cos(theta));
733 return beta;
734}
735
736
737
738
739void AliRICHPatRec::HoughResponse()
740
741{
237c933d 742
743// Implement Hough response pat. rec. method
744
e7257cad 745 int bin=0;
746 int bin1=0;
747 int bin2=0;
237c933d 748 int i, j, k, nCorrBand;
749 int etaBin = 750;
750 float hcs[750];
751 float angle, thetaCerMean;
e7257cad 752
237c933d 753 float etaPeak[30];
754 float etaMin = 0.00;
755 float etaMax = 0.75;
756 float stepEta = 0.001;
757 float windowEta = 0.040;
e7257cad 758
237c933d 759 int nBin;
e7257cad 760
237c933d 761 float etaPeakPos = -1;
762 Int_t etaPeakCount = -1;
e7257cad 763
237c933d 764 thetaCerMean = 0.;
e7257cad 765 fThetaCerenkov = 0.;
766
237c933d 767 nBin = (int)(0.5+etaMax/(stepEta));
768 nCorrBand = (int)(0.5+ windowEta/(2 * stepEta));
769 memset ((void *)hcs, 0, etaBin*sizeof(int));
e7257cad 770
771 for (k=0; k< fNumEtaPhotons; k++) {
772
773 angle = fEtaPhotons[k];
774
237c933d 775 if (angle>=etaMin && angle<= etaMax) {
776 bin = (int)(0.5+angle/(stepEta));
777 bin1= bin-nCorrBand;
778 bin2= bin+nCorrBand;
e7257cad 779 if (bin1<0) bin1=0;
237c933d 780 if (bin2>nBin) bin2=nBin;
e7257cad 781
782 for (j=bin1; j<bin2; j++) {
237c933d 783 hcs[j] += fWeightPhotons[k];
e7257cad 784 }
785
237c933d 786 thetaCerMean += angle;
e7257cad 787 }
788 }
789
237c933d 790 thetaCerMean /= fNumEtaPhotons;
e7257cad 791
237c933d 792 HoughFiltering(hcs);
793
794 for (bin=0; bin <nBin; bin++) {
795 angle = (bin+0.5) * (stepEta);
796 if (hcs[bin] && hcs[bin] > etaPeakPos) {
797 etaPeakCount = 0;
798 etaPeakPos = hcs[bin];
799 etaPeak[0]=angle;
e7257cad 800 }
801 else {
237c933d 802 if (hcs[bin] == etaPeakPos) {
803 etaPeak[++etaPeakCount] = angle;
e7257cad 804 }
805 }
806 }
807
237c933d 808 for (i=0; i<etaPeakCount+1; i++) {
809 fThetaCerenkov += etaPeak[i];
e7257cad 810 }
237c933d 811 if (etaPeakCount>=0) {
812 fThetaCerenkov /= etaPeakCount+1;
813 fThetaPeakPos = etaPeakPos;
e7257cad 814 }
815}
816
817
237c933d 818void AliRICHPatRec::HoughFiltering(float hcs[])
e7257cad 819{
237c933d 820
821// hough filtering
822
823 float hcsFilt[750];
824 float k[5] = {0.05, 0.25, 0.4, 0.25, 0.05};
825 int nx, i, nxDx;
e7257cad 826 int sizeHCS;
237c933d 827 int nBin;
e7257cad 828
237c933d 829 int etaBin = 750;
830 float etaMax = 0.75;
831 float stepEta = 0.001;
e7257cad 832
237c933d 833 nBin = (int)(1+etaMax/stepEta);
834 sizeHCS = etaBin*sizeof(float);
e7257cad 835
237c933d 836 memset ((void *)hcsFilt, 0, sizeHCS);
e7257cad 837
237c933d 838 for (nx = 0; nx < nBin; nx++) {
e7257cad 839 for (i = 0; i < 5; i++) {
237c933d 840 nxDx = nx + (i-2);
841 if (nxDx> -1 && nxDx<nBin)
842 hcsFilt[nx] += hcs[nxDx] * k[i];
e7257cad 843 }
844 }
845
237c933d 846 for (nx = 0; nx < nBin; nx++) {
847 hcs[nx] = hcsFilt[nx];
e7257cad 848 }
849}
850
e0b63a71 851/*void AliRICHPatRec::CerenkovRingDrawing()
e7257cad 852{
853
854//to draw Cherenkov ring by known Cherenkov angle
855
e0b63a71 856 Int_t nmaxdegrees;
857 Int_t Nphpad;
858 Float_t phpad;
237c933d 859 Float_t nfreonave, nquartzave;
860 Float_t aveEnerg;
861 Float_t energy[2];
862 Float_t e1, e2, f1, f2;
863 Float_t bandradius;
e0b63a71 864
e7257cad 865//parameters to calculate freon window refractive index vs. energy
e0b63a71 866
237c933d 867 Float_t a = 1.177;
868 Float_t b = 0.0172;
869
e7257cad 870//parameters to calculate quartz window refractive index vs. energy
e0b63a71 871
237c933d 872 energy[0] = 5.0;
873 energy[1] = 8.0;
874 e1 = 10.666;
875 e2 = 18.125;
876 f1 = 46.411;
877 f2 = 228.71;
878
e7257cad 879
e0b63a71 880 nmaxdegrees = 36;
237c933d 881
e0b63a71 882 for (Nphpad=0; Nphpad<nmaxdegrees;Nphpad++) {
883
884 phpad = (360./(Float_t)nmaxdegrees)*(Float_t)Nphpad;
e7257cad 885
237c933d 886 aveEnerg = (energy[0]+energy[1])/2.;
887
888 nfreonave = a+b*aveEnerg;
00df6e79 889 nquartzave = sqrt(1+(f1/(TMath::Power(e1,2)-TMath::Power(aveEnerg,2)))+
890 (f2/(TMath::Power(e2,2)-TMath::Power(aveEnerg,2))));
237c933d 891
237c933d 892 bandradius = DistanceFromMip(nfreonave, nquartzave,
e0b63a71 893 fEmissPoint,fThetaCerenkov, phpad);
e7257cad 894
e0b63a71 895 fCoordEllipse[0][Nphpad] = fOnCathode[0];
896 fCoordEllipse[1][Nphpad] = fOnCathode[1];
897 printf(" values %f %f \n",fOnCathode[0],fOnCathode[1]);
237c933d 898
237c933d 899 }
900
e0b63a71 901}*/
e7257cad 902
e7257cad 903