]>
Commit | Line | Data |
---|---|---|
51ad6848 | 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 | /* $Id$ */ | |
17 | ||
18 | //------------------------------------------------------------------------- | |
c1e38247 | 19 | // |
20 | // Implementation of the ESD V0MI vertex class | |
21 | // This class is part of the Event Data Summary | |
22 | // set of classes and contains information about | |
23 | // V0 kind vertexes generated by a neutral particle | |
24 | // Numerical part - AliHelix functionality used | |
25 | // | |
26 | // Likelihoods for Point angle, DCA and Causality defined => can be used as cut parameters | |
27 | // HIGHLY recomended | |
28 | // | |
29 | // Quality information can be used as additional cut variables | |
30 | // | |
51ad6848 | 31 | // Origin: Marian Ivanov marian.ivanov@cern.ch |
32 | //------------------------------------------------------------------------- | |
33 | ||
34 | #include <Riostream.h> | |
35 | #include <TMath.h> | |
0703142d | 36 | |
51ad6848 | 37 | #include "AliESDV0MI.h" |
38 | #include "AliHelix.h" | |
39 | ||
40 | ||
41 | ClassImp(AliESDV0MI) | |
42 | ||
c1e38247 | 43 | AliESDV0MIParams AliESDV0MI::fgkParams; |
44 | ||
45 | ||
90e48c0c | 46 | AliESDV0MI::AliESDV0MI() : |
47 | AliESDv0(), | |
48 | fParamP(), | |
49 | fParamM(), | |
50 | fID(0), | |
51 | fDist1(-1), | |
52 | fDist2(-1), | |
53 | fRr(-1), | |
54 | fStatus(0), | |
55 | fRow0(-1), | |
56 | fDistNorm(0), | |
57 | fDistSigma(0), | |
58 | fChi2Before(0), | |
59 | fNBefore(0), | |
60 | fChi2After(0), | |
61 | fNAfter(0), | |
62 | fPointAngleFi(0), | |
63 | fPointAngleTh(0), | |
64 | fPointAngle(0) | |
65 | { | |
51ad6848 | 66 | // |
67 | //Dafault constructor | |
68 | // | |
81e97e0d | 69 | for (Int_t i=0;i<4;i++){fCausality[i]=0;} |
6605de26 | 70 | for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;} |
29641977 | 71 | for (Int_t i=0;i<2;i++){fNormDCAPrim[0]=0;fNormDCAPrim[1]=0;} |
81e97e0d | 72 | } |
73 | ||
c1e38247 | 74 | Double_t AliESDV0MI::GetSigmaY(){ |
75 | // | |
76 | // return sigmay in y at vertex position using covariance matrix | |
77 | // | |
78 | const Double_t * cp = fParamP.GetCovariance(); | |
79 | const Double_t * cm = fParamM.GetCovariance(); | |
80 | Double_t sigmay = cp[0]+cm[0]+ cp[5]*(fParamP.X()-fRr)*(fParamP.X()-fRr)+ cm[5]*(fParamM.X()-fRr)*(fParamM.X()-fRr); | |
81 | return (sigmay>0) ? TMath::Sqrt(sigmay):100; | |
82 | } | |
83 | ||
84 | Double_t AliESDV0MI::GetSigmaZ(){ | |
85 | // | |
86 | // return sigmay in y at vertex position using covariance matrix | |
87 | // | |
88 | const Double_t * cp = fParamP.GetCovariance(); | |
89 | const Double_t * cm = fParamM.GetCovariance(); | |
90 | Double_t sigmaz = cp[2]+cm[2]+ cp[9]*(fParamP.X()-fRr)*(fParamP.X()-fRr)+ cm[9]*(fParamM.X()-fRr)*(fParamM.X()-fRr); | |
91 | return (sigmaz>0) ? TMath::Sqrt(sigmaz):100; | |
92 | } | |
93 | ||
94 | Double_t AliESDV0MI::GetSigmaD0(){ | |
95 | // | |
96 | // Sigma parameterization using covariance matrix | |
97 | // | |
98 | // sigma of distance between two tracks in vertex position | |
99 | // sigma of DCA is proportianal to sigmaD0 | |
100 | // factor 2 difference is explained by the fact that the DCA is calculated at the position | |
101 | // where the tracks as closest together ( not exact position of the vertex) | |
102 | // | |
103 | const Double_t * cp = fParamP.GetCovariance(); | |
104 | const Double_t * cm = fParamM.GetCovariance(); | |
105 | Double_t sigmaD0 = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0; | |
106 | sigmaD0 += ((fParamP.X()-fRr)*(fParamP.X()-fRr))*(cp[5]+cp[9]); | |
107 | sigmaD0 += ((fParamM.X()-fRr)*(fParamM.X()-fRr))*(cm[5]+cm[9]); | |
108 | return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100; | |
109 | } | |
110 | ||
111 | ||
112 | Double_t AliESDV0MI::GetSigmaAP0(){ | |
113 | // | |
114 | //Sigma parameterization using covariance matrices | |
115 | // | |
116 | Double_t prec = TMath::Sqrt((fPM[0]+fPP[0])*(fPM[0]+fPP[0]) | |
117 | +(fPM[1]+fPP[1])*(fPM[1]+fPP[1]) | |
118 | +(fPM[2]+fPP[2])*(fPM[2]+fPP[2])); | |
119 | Double_t normp = TMath::Sqrt(fPP[0]*fPP[0]+fPP[1]*fPP[1]+fPP[2]*fPP[2])/prec; // fraction of the momenta | |
120 | Double_t normm = TMath::Sqrt(fPM[0]*fPM[0]+fPM[1]*fPM[1]+fPM[2]*fPM[2])/prec; | |
121 | const Double_t * cp = fParamP.GetCovariance(); | |
122 | const Double_t * cm = fParamM.GetCovariance(); | |
123 | Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0; // minimal part | |
124 | sigmaAP0 += (cp[5]+cp[9])*(normp*normp)+(cm[5]+cm[9])*(normm*normm); // angular resolution part | |
125 | Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01); // vertex position part | |
126 | sigmaAP0 += 0.5*sigmaAP1*sigmaAP1; | |
127 | return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100; | |
128 | } | |
129 | ||
130 | Double_t AliESDV0MI::GetEffectiveSigmaD0(){ | |
131 | // | |
132 | // minimax - effective Sigma parameterization | |
133 | // p12 effective curvature and v0 radius postion used as parameters | |
134 | // | |
135 | Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+ | |
136 | fParamM.GetParameter()[4]*fParamM.GetParameter()[4]); | |
137 | Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)*fgkParams.fPSigmaCoefDE*p12*p12; | |
138 | sigmaED0*= sigmaED0; | |
139 | sigmaED0*= sigmaED0; | |
140 | sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE); | |
141 | return (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE; | |
142 | } | |
143 | ||
144 | ||
145 | Double_t AliESDV0MI::GetEffectiveSigmaAP0(){ | |
146 | // | |
147 | // effective Sigma parameterization of point angle resolution | |
148 | // | |
149 | Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+ | |
150 | fParamM.GetParameter()[4]*fParamM.GetParameter()[4]); | |
151 | Double_t sigmaAPE= fgkParams.fPSigmaBase0APE; | |
152 | sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr); | |
153 | sigmaAPE*= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APE*p12); | |
154 | sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE); | |
155 | return sigmaAPE; | |
156 | } | |
157 | ||
158 | ||
159 | Double_t AliESDV0MI::GetMinimaxSigmaAP0(){ | |
160 | // | |
161 | // calculate mini-max effective sigma of point angle resolution | |
162 | // | |
163 | //compv0->fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)*0.5,1.5*SigmaAPE0),0.5*SigmaAPE0+0.003)"); | |
164 | Double_t effectiveSigma = GetEffectiveSigmaAP0(); | |
165 | Double_t sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma); | |
166 | sigmaMMAP = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma); | |
167 | sigmaMMAP = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0); | |
168 | return sigmaMMAP; | |
169 | } | |
170 | Double_t AliESDV0MI::GetMinimaxSigmaD0(){ | |
171 | // | |
172 | // calculate mini-max sigma of dca resolution | |
173 | // | |
174 | //compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)*0.5,1.5*SigmaDE0),0.5*SigmaDE0)"); | |
175 | Double_t effectiveSigma = GetEffectiveSigmaD0(); | |
176 | Double_t sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma); | |
177 | sigmaMMD0 = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma); | |
178 | sigmaMMD0 = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0); | |
179 | return sigmaMMD0; | |
180 | } | |
181 | ||
182 | ||
183 | Double_t AliESDV0MI::GetLikelihoodAP(Int_t mode0, Int_t mode1){ | |
184 | // | |
185 | // get likelihood for point angle | |
186 | // | |
187 | Double_t sigmaAP = 0.007; //default sigma | |
188 | switch (mode0){ | |
189 | case 0: | |
190 | sigmaAP = GetSigmaAP0(); // mode 0 - covariance matrix estimates used | |
191 | break; | |
192 | case 1: | |
193 | sigmaAP = GetEffectiveSigmaAP0(); // mode 1 - effective sigma used | |
194 | break; | |
195 | case 2: | |
196 | sigmaAP = GetMinimaxSigmaAP0(); // mode 2 - minimax sigma | |
197 | break; | |
198 | } | |
199 | Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.); | |
200 | //normalized point angle, restricted - because of overflow problems in Exp | |
201 | Double_t likelihood = 0; | |
202 | switch(mode1){ | |
203 | case 0: | |
204 | likelihood = TMath::Exp(-0.5*apNorm*apNorm); | |
205 | // one component | |
206 | break; | |
207 | case 1: | |
208 | likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm))/1.5; | |
209 | // two components | |
210 | break; | |
211 | case 2: | |
212 | likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm)+0.25*TMath::Exp(-0.125*apNorm*apNorm))/1.75; | |
213 | // three components | |
214 | break; | |
215 | } | |
216 | return likelihood; | |
217 | } | |
218 | ||
219 | Double_t AliESDV0MI::GetLikelihoodD(Int_t mode0, Int_t mode1){ | |
220 | // | |
221 | // get likelihood for DCA | |
222 | // | |
223 | Double_t sigmaD = 0.03; //default sigma | |
224 | switch (mode0){ | |
225 | case 0: | |
226 | sigmaD = GetSigmaD0(); // mode 0 - covariance matrix estimates used | |
227 | break; | |
228 | case 1: | |
229 | sigmaD = GetEffectiveSigmaD0(); // mode 1 - effective sigma used | |
230 | break; | |
231 | case 2: | |
232 | sigmaD = GetMinimaxSigmaD0(); // mode 2 - minimax sigma | |
233 | break; | |
234 | } | |
235 | Double_t dNorm = TMath::Min(fDist2/sigmaD,50.); | |
236 | //normalized point angle, restricted - because of overflow problems in Exp | |
237 | Double_t likelihood = 0; | |
238 | switch(mode1){ | |
239 | case 0: | |
240 | likelihood = TMath::Exp(-2.*dNorm); | |
241 | // one component | |
242 | break; | |
243 | case 1: | |
244 | likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm))/1.5; | |
245 | // two components | |
246 | break; | |
247 | case 2: | |
248 | likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm)+0.25*TMath::Exp(-0.5*dNorm))/1.75; | |
249 | // three components | |
250 | break; | |
251 | } | |
252 | return likelihood; | |
253 | ||
254 | } | |
255 | ||
256 | Double_t AliESDV0MI::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/){ | |
257 | // | |
258 | // get likelihood for Causality | |
259 | // !!! Causality variables defined in AliITStrackerMI !!! | |
260 | // when more information was available | |
261 | // | |
262 | Double_t likelihood = 0.5; | |
263 | Double_t minCausal = TMath::Min(fCausality[0],fCausality[1]); | |
264 | Double_t maxCausal = TMath::Max(fCausality[0],fCausality[1]); | |
265 | // minCausal = TMath::Max(minCausal,0.5*maxCausal); | |
266 | //compv0->fTree->SetAlias("LCausal","(1.05-(2*(0.8-exp(-max(RC.fV0rec.fCausality[0],RC.fV0rec.fCausality[1])))+2*(0.8-exp(-min(RC.fV0rec.fCausality[0],RC.fV0rec.fCausality[1]))))/2)**4"); | |
267 | ||
268 | switch(mode0){ | |
269 | case 0: | |
270 | //normalization | |
271 | likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.); | |
272 | break; | |
273 | case 1: | |
274 | likelihood = TMath::Power(1.05-(2*(0.8-TMath::Exp(-maxCausal))+(2*(0.8-TMath::Exp(-minCausal))))*0.5,4.); | |
275 | break; | |
276 | } | |
277 | return likelihood; | |
278 | ||
279 | } | |
81e97e0d | 280 | |
281 | void AliESDV0MI::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1) | |
282 | { | |
283 | // | |
284 | // set probabilities | |
285 | // | |
286 | fCausality[0] = pb0; // probability - track 0 exist before vertex | |
287 | fCausality[1] = pb1; // probability - track 1 exist before vertex | |
288 | fCausality[2] = pa0; // probability - track 0 exist close after vertex | |
289 | fCausality[3] = pa1; // probability - track 1 exist close after vertex | |
51ad6848 | 290 | } |
6605de26 | 291 | void AliESDV0MI::SetClusters(Int_t *clp, Int_t *clm) |
292 | { | |
293 | // | |
294 | // Set its clusters indexes | |
295 | // | |
296 | for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i]; | |
297 | for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i]; | |
298 | } | |
299 | ||
51ad6848 | 300 | |
301 | void AliESDV0MI::SetP(const AliExternalTrackParam & paramp) { | |
302 | // | |
81e97e0d | 303 | // set track + |
51ad6848 | 304 | // |
305 | fParamP = paramp; | |
306 | } | |
307 | ||
308 | void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){ | |
309 | // | |
81e97e0d | 310 | //set track - |
51ad6848 | 311 | // |
312 | fParamM = paramm; | |
51ad6848 | 313 | } |
314 | ||
81e97e0d | 315 | void AliESDV0MI::SetRp(const Double_t *rp){ |
316 | // | |
317 | // set pid + | |
318 | // | |
319 | for (Int_t i=0;i<5;i++) fRP[i]=rp[i]; | |
320 | } | |
321 | ||
322 | void AliESDV0MI::SetRm(const Double_t *rm){ | |
323 | // | |
324 | // set pid - | |
325 | // | |
326 | for (Int_t i=0;i<5;i++) fRM[i]=rm[i]; | |
327 | } | |
328 | ||
329 | ||
51ad6848 | 330 | void AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5]) |
331 | { | |
332 | // | |
333 | // set PID hypothesy | |
334 | // | |
335 | // norm PID to 1 | |
336 | Float_t sump =0; | |
337 | Float_t summ =0; | |
338 | for (Int_t i=0;i<5;i++){ | |
339 | fRP[i]=pidp[i]; | |
340 | sump+=fRP[i]; | |
341 | fRM[i]=pidm[i]; | |
342 | summ+=fRM[i]; | |
343 | } | |
344 | for (Int_t i=0;i<5;i++){ | |
345 | fRP[i]/=sump; | |
346 | fRM[i]/=summ; | |
347 | } | |
348 | } | |
349 | ||
350 | Float_t AliESDV0MI::GetProb(UInt_t p1, UInt_t p2){ | |
351 | // | |
352 | // | |
353 | // | |
354 | // | |
355 | return TMath::Max(fRP[p1]+fRM[p2], fRP[p2]+fRM[p1]); | |
356 | } | |
357 | ||
358 | Float_t AliESDV0MI::GetEffMass(UInt_t p1, UInt_t p2){ | |
359 | // | |
360 | // calculate effective mass | |
361 | // | |
0703142d | 362 | const Float_t kpmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01, |
51ad6848 | 363 | 4.93599999999999983e-01, 9.38270000000000048e-01}; |
364 | if (p1>4) return -1; | |
365 | if (p2>4) return -1; | |
0703142d | 366 | Float_t mass1 = kpmass[p1]; |
367 | Float_t mass2 = kpmass[p2]; | |
51ad6848 | 368 | Double_t *m1 = fPP; |
369 | Double_t *m2 = fPM; | |
370 | // | |
6605de26 | 371 | //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){ |
372 | // m1 = fPM; | |
373 | // m2 = fPP; | |
374 | //} | |
51ad6848 | 375 | // |
376 | Float_t e1 = TMath::Sqrt(mass1*mass1+ | |
377 | m1[0]*m1[0]+ | |
378 | m1[1]*m1[1]+ | |
379 | m1[2]*m1[2]); | |
380 | Float_t e2 = TMath::Sqrt(mass2*mass2+ | |
381 | m2[0]*m2[0]+ | |
382 | m2[1]*m2[1]+ | |
383 | m2[2]*m2[2]); | |
384 | Float_t mass = | |
385 | (m2[0]+m1[0])*(m2[0]+m1[0])+ | |
386 | (m2[1]+m1[1])*(m2[1]+m1[1])+ | |
387 | (m2[2]+m1[2])*(m2[2]+m1[2]); | |
388 | ||
389 | mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass); | |
390 | return mass; | |
391 | } | |
392 | ||
393 | void AliESDV0MI::Update(Float_t vertex[3]) | |
394 | { | |
395 | // | |
396 | // updates Kink Info | |
397 | // | |
81e97e0d | 398 | // Float_t distance1,distance2; |
399 | Float_t distance2; | |
51ad6848 | 400 | // |
401 | AliHelix phelix(fParamP); | |
402 | AliHelix mhelix(fParamM); | |
403 | // | |
404 | //find intersection linear | |
405 | // | |
406 | Double_t phase[2][2],radius[2]; | |
407 | Int_t points = phelix.GetRPHIintersections(mhelix, phase, radius,200); | |
408 | Double_t delta1=10000,delta2=10000; | |
81e97e0d | 409 | /* |
b07a036b | 410 | if (points<=0) return; |
51ad6848 | 411 | if (points>0){ |
412 | phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
413 | phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
414 | phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
415 | } | |
416 | if (points==2){ | |
417 | phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
418 | phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
419 | phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
420 | } | |
421 | distance1 = TMath::Min(delta1,delta2); | |
81e97e0d | 422 | */ |
51ad6848 | 423 | // |
424 | //find intersection parabolic | |
425 | // | |
426 | points = phelix.GetRPHIintersections(mhelix, phase, radius); | |
427 | delta1=10000,delta2=10000; | |
428 | Double_t d1=1000.,d2=10000.; | |
29641977 | 429 | Double_t err[3],angles[3]; |
b07a036b | 430 | if (points<=0) return; |
51ad6848 | 431 | if (points>0){ |
432 | phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
433 | phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); | |
29641977 | 434 | if (TMath::Abs(fParamP.X()-TMath::Sqrt(radius[0])<3) && TMath::Abs(fParamM.X()-TMath::Sqrt(radius[0])<3)){ |
435 | // if we are close to vertex use error parama | |
436 | // | |
437 | err[1] = fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]+0.05*0.05 | |
438 | +0.3*(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]); | |
439 | err[2] = fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]+0.05*0.05 | |
440 | +0.3*(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]); | |
441 | ||
442 | phelix.GetAngle(phase[0][0],mhelix,phase[0][1],angles); | |
443 | Double_t tfi = TMath::Abs(TMath::Tan(angles[0])); | |
444 | Double_t tlam = TMath::Abs(TMath::Tan(angles[1])); | |
445 | err[0] = err[1]/((0.2+tfi)*(0.2+tfi))+err[2]/((0.2+tlam)*(0.2+tlam)); | |
446 | err[0] = ((err[1]*err[2]/((0.2+tfi)*(0.2+tfi)*(0.2+tlam)*(0.2+tlam))))/err[0]; | |
447 | phelix.ParabolicDCA2(mhelix,phase[0][0],phase[0][1],radius[0],delta1,err); | |
448 | } | |
51ad6848 | 449 | Double_t xd[3],xm[3]; |
450 | phelix.Evaluate(phase[0][0],xd); | |
451 | mhelix.Evaluate(phase[0][1],xm); | |
452 | d1 = (xd[0]-xm[0])*(xd[0]-xm[0])+(xd[1]-xm[1])*(xd[1]-xm[1])+(xd[2]-xm[2])*(xd[2]-xm[2]); | |
453 | } | |
454 | if (points==2){ | |
455 | phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
456 | phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); | |
29641977 | 457 | if (TMath::Abs(fParamP.X()-TMath::Sqrt(radius[1])<3) && TMath::Abs(fParamM.X()-TMath::Sqrt(radius[1])<3)){ |
458 | // if we are close to vertex use error paramatrization | |
459 | // | |
460 | err[1] = fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]+0.05*0.05 | |
461 | +0.3*(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]); | |
462 | err[2] = fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]+0.05*0.05 | |
463 | +0.3*(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]); | |
464 | ||
465 | phelix.GetAngle(phase[1][0],mhelix,phase[1][1],angles); | |
466 | Double_t tfi = TMath::Abs(TMath::Tan(angles[0])); | |
467 | Double_t tlam = TMath::Abs(TMath::Tan(angles[1])); | |
468 | err[0] = err[1]/((0.2+tfi)*(0.2+tfi))+err[2]/((0.2+tlam)*(0.2+tlam)); | |
469 | err[0] = ((err[1]*err[2]/((0.2+tfi)*(0.2+tfi)*(0.2+tlam)*(0.2+tlam))))/err[0]; | |
470 | phelix.ParabolicDCA2(mhelix,phase[1][0],phase[1][1],radius[1],delta2,err); | |
471 | } | |
51ad6848 | 472 | Double_t xd[3],xm[3]; |
473 | phelix.Evaluate(phase[1][0],xd); | |
474 | mhelix.Evaluate(phase[1][1],xm); | |
475 | d2 = (xd[0]-xm[0])*(xd[0]-xm[0])+(xd[1]-xm[1])*(xd[1]-xm[1])+(xd[2]-xm[2])*(xd[2]-xm[2]); | |
476 | } | |
477 | // | |
478 | distance2 = TMath::Min(delta1,delta2); | |
479 | if (delta1<delta2){ | |
480 | //get V0 info | |
481 | Double_t xd[3],xm[3]; | |
482 | phelix.Evaluate(phase[0][0],xd); | |
483 | mhelix.Evaluate(phase[0][1], xm); | |
484 | fXr[0] = 0.5*(xd[0]+xm[0]); | |
485 | fXr[1] = 0.5*(xd[1]+xm[1]); | |
486 | fXr[2] = 0.5*(xd[2]+xm[2]); | |
29641977 | 487 | |
488 | Float_t wy = fParamP.GetCovariance()[0]/(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]); | |
489 | Float_t wz = fParamP.GetCovariance()[2]/(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]); | |
490 | fXr[0] = 0.5*( (1.-wy)*xd[0]+ wy*xm[0] + (1.-wz)*xd[0]+ wz*xm[0] ); | |
491 | fXr[1] = (1.-wy)*xd[1]+ wy*xm[1]; | |
492 | fXr[2] = (1.-wz)*xd[2]+ wz*xm[2]; | |
51ad6848 | 493 | // |
494 | phelix.GetMomentum(phase[0][0],fPP); | |
495 | mhelix.GetMomentum(phase[0][1],fPM); | |
496 | phelix.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle); | |
497 | fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]); | |
498 | } | |
499 | else{ | |
500 | Double_t xd[3],xm[3]; | |
501 | phelix.Evaluate(phase[1][0],xd); | |
502 | mhelix.Evaluate(phase[1][1], xm); | |
503 | fXr[0] = 0.5*(xd[0]+xm[0]); | |
504 | fXr[1] = 0.5*(xd[1]+xm[1]); | |
505 | fXr[2] = 0.5*(xd[2]+xm[2]); | |
29641977 | 506 | Float_t wy = fParamP.GetCovariance()[0]/(fParamP.GetCovariance()[0]+fParamM.GetCovariance()[0]); |
507 | Float_t wz = fParamP.GetCovariance()[2]/(fParamP.GetCovariance()[2]+fParamM.GetCovariance()[2]); | |
508 | fXr[0] = 0.5*( (1.-wy)*xd[0]+ wy*xm[0] + (1.-wz)*xd[0]+ wz*xm[0] ); | |
509 | fXr[1] = (1.-wy)*xd[1]+ wy*xm[1]; | |
510 | fXr[2] = (1.-wz)*xd[2]+ wz*xm[2]; | |
51ad6848 | 511 | // |
512 | phelix.GetMomentum(phase[1][0], fPP); | |
513 | mhelix.GetMomentum(phase[1][1], fPM); | |
514 | phelix.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle); | |
515 | fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]); | |
516 | } | |
517 | fDist1 = TMath::Sqrt(TMath::Min(d1,d2)); | |
518 | fDist2 = TMath::Sqrt(distance2); | |
519 | // | |
520 | // | |
81e97e0d | 521 | Double_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]}; |
522 | Double_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]}; | |
523 | Double_t vnorm2 = v[0]*v[0]+v[1]*v[1]; | |
c1e38247 | 524 | if (TMath::Abs(v[2])>100000) return; |
525 | Double_t vnorm3 = TMath::Sqrt(TMath::Abs(v[2]*v[2]+vnorm2)); | |
51ad6848 | 526 | vnorm2 = TMath::Sqrt(vnorm2); |
81e97e0d | 527 | Double_t pnorm2 = p[0]*p[0]+p[1]*p[1]; |
528 | Double_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2); | |
51ad6848 | 529 | pnorm2 = TMath::Sqrt(pnorm2); |
530 | fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2); | |
531 | fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3); | |
532 | fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3); | |
533 | // | |
534 | } | |
535 |