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