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