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