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