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