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