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