]>
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$ | |
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 | ||
88 | ClassImp(AliRICHPatRec) | |
89 | //___________________________________________ | |
90 | AliRICHPatRec::AliRICHPatRec() : TObject() | |
91 | { | |
237c933d | 92 | // Default constructor |
93 | ||
e7257cad | 94 | //fChambers = 0; |
95 | } | |
96 | //___________________________________________ | |
97 | AliRICHPatRec::AliRICHPatRec(const char *name, const char *title) | |
98 | : TObject() | |
99 | { | |
237c933d | 100 | //Constructor for Bari's pattern recogniton method object |
e7257cad | 101 | } |
102 | ||
103 | void 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 | 211 | Int_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 | ||
288 | Float_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 | ||
313 | Float_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 | ||
381 | void 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 | ||
398 | void 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 | ||
406 | void 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 | ||
516 | void 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 | ||
531 | Int_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 | ||
600 | Float_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 | 698 | Float_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 | ||
732 | Float_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 | ||
751 | Float_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 | ||
767 | Float_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 | ||
781 | void 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 | 860 | void 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 |