]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliRecInfo.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
inconsistency corrected, version introduced
[u/mrichter/AliRoot.git] / PWG1 / AliRecInfo.cxx
CommitLineData
06d06fbb 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///////////////////////////////////////////////////////////////////////////////
18// //
19// Time Projection Chamber //
20// Comparison macro for reconstructed tracks - ESDs V0s //
21// responsible:
22// marian.ivanov@cern.ch //
23//
24//
25
26
27
28
29
30#include <stdio.h>
31#include <string.h>
32//ROOT includes
33#include "Rtypes.h"
34#include "TFile.h"
35#include "TTree.h"
36#include "TChain.h"
37#include "TCut.h"
38#include "TString.h"
39#include "TBenchmark.h"
40#include "TStopwatch.h"
41#include "TParticle.h"
42#include "TSystem.h"
43#include "TTimer.h"
44#include "TVector3.h"
45#include "TPad.h"
46#include "TCanvas.h"
47#include "TH1F.h"
48#include "TH2F.h"
49#include "TF1.h"
50#include "TText.h"
51#include "Getline.h"
52#include "TStyle.h"
53//ALIROOT includes
54#include "AliRun.h"
55#include "AliStack.h"
56#include "AliESDtrack.h"
57#include "AliSimDigits.h"
58#include "AliTPCParam.h"
59#include "AliTPC.h"
60#include "AliTPCLoader.h"
61#include "AliDetector.h"
62#include "AliTrackReference.h"
63#include "AliRun.h"
64#include "AliTPCParamSR.h"
65#include "AliTracker.h"
66#include "AliComplexCluster.h"
67#include "AliMagF.h"
68#include "AliESD.h"
69#include "AliESDfriend.h"
70#include "AliESDtrack.h"
71#include "AliTPCseed.h"
72#include "AliITStrackMI.h"
73#include "AliTRDtrack.h"
74#include "AliHelix.h"
75#include "AliESDVertex.h"
76#include "AliExternalTrackParam.h"
77#include "AliESDkink.h"
78#include "AliESDv0.h"
79#include "AliV0.h"
80//
81#include "AliTreeDraw.h"
82#include "AliGenInfo.h"
83#include "AliRecInfo.h"
84
85
86
87ClassImp(AliESDRecInfo)
88ClassImp(AliESDRecV0Info)
89ClassImp(AliESDRecKinkInfo)
90
91
92
93
94AliTPCParam * GetTPCParam(){
95 AliTPCParamSR * par = new AliTPCParamSR;
96 par->Update();
97 return par;
98}
99
100
101
102
103AliESDRecInfo::AliESDRecInfo():
104 fITSOn(0), // ITS refitted inward
105 fTRDOn(0), // ITS refitted inward
106 fDeltaP(0), //delta of momenta
107 fSign(0), // sign
108 fReconstructed(0), //flag if track was reconstructed
109 fFake(0), // fake track
110 fMultiple(0), // number of reconstructions
111 fTPCOn(0), // TPC refitted inward
112 fBestTOFmatch(0), //best matching between times
113 fESDtrack(0), // esd track
114 fTrackF(0), // friend track
115 fTPCtrack(0), // tpc track
116 fITStrack(0), // its track
117 fTRDtrack(0) // trd track
118{
119 //
120 // default constructor
121 //
122}
123
124
125AliESDRecInfo::AliESDRecInfo(const AliESDRecInfo& recinfo):
126 TObject()
127{
128 //
129 //
130 //
131 memcpy(this,&recinfo, sizeof(recinfo));
132 fESDtrack=0; fTrackF=0; fTPCtrack=0;fITStrack=0;fTRDtrack=0;
133 SetESDtrack(recinfo.GetESDtrack());
134}
135
136
137AliESDRecInfo::~AliESDRecInfo()
138
139{
140 //
141 // destructor
142 //
143 if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
144 if (fTrackF) { delete fTrackF; fTrackF=0;}
145 if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
146 if (fITStrack) { delete fITStrack; fITStrack=0;}
147 if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
148
149}
150
151
152
153void AliESDRecInfo::Reset()
154{
155 //
156 // reset info
157 //
158 fMultiple =0;
159 fFake =0;
160 fReconstructed=0;
161 if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
162 if (fTrackF) { delete fTrackF; fTrackF=0;}
163 if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
164 if (fITStrack) { delete fITStrack; fITStrack=0;}
165 if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
166}
167
168void AliESDRecInfo::SetESDtrack(const AliESDtrack *track){
169 //
170 //
171 //
172 if (fESDtrack) delete fESDtrack;
173 fESDtrack = (AliESDtrack*)track->Clone();
174 if (track->GetFriendTrack()){
175 if (fTrackF) delete fTrackF;
176 fTrackF = (AliESDfriendTrack*)track->GetFriendTrack()->Clone();
177 if (fTrackF->GetCalibObject(0)){
178 if (fTPCtrack) delete fTPCtrack;
179 fTPCtrack = (AliTPCseed*)fTrackF->GetCalibObject(0)->Clone();
180 }
181 }
182
183}
184
185void AliESDRecInfo::UpdatePoints(AliESDtrack*track)
186{
187 //
188 //
189 Int_t iclusters[200];
190 Float_t density[160];
191 for (Int_t i=0;i<160;i++) density[i]=-1.;
192 fTPCPoints[0]= 160;
193 fTPCPoints[1] = -1;
194 //
195 if (fTPCPoints[0]<fTPCPoints[1]) return;
196 // Int_t nclusters=track->GetTPCclusters(iclusters);
197
198 Int_t ngood=0;
199 Int_t undeff=0;
200 Int_t nall =0;
201 Int_t range=20;
202 for (Int_t i=0;i<160;i++){
203 Int_t last = i-range;
204 if (nall<range) nall++;
205 if (last>=0){
206 if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--;
207 if (iclusters[last]==-1) undeff--;
208 }
209 if (iclusters[i]>0&& (iclusters[i]&0x8000)==0) ngood++;
210 if (iclusters[i]==-1) undeff++;
211 if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
212 }
213 Float_t maxdens=0;
214 Int_t indexmax =0;
215 for (Int_t i=0;i<160;i++){
216 if (density[i]<0) continue;
217 if (density[i]>maxdens){
218 maxdens=density[i];
219 indexmax=i;
220 }
221 }
222 //
223 //max dens point
224 fTPCPoints[3] = maxdens;
225 fTPCPoints[1] = indexmax;
226 //
227 // last point
228 for (Int_t i=indexmax;i<160;i++){
229 if (density[i]<0) continue;
230 if (density[i]<maxdens/2.) {
231 break;
232 }
233 fTPCPoints[2]=i;
234 }
235 //
236 // first point
237 for (Int_t i=indexmax;i>0;i--){
238 if (density[i]<0) continue;
239 if (density[i]<maxdens/2.) {
240 break;
241 }
242 fTPCPoints[0]=i;
243 }
244 //
245 // Density at the last 30 padrows
246 //
247 //
248 nall = 0;
249 ngood = 0;
250 for (Int_t i=159;i>0;i--){
251 if (iclusters[i]==-1) continue; //dead zone
252 nall++;
253 if (iclusters[i]>0) ngood++;
254 if (nall>20) break;
255 }
256 fTPCPoints[4] = Float_t(ngood)/Float_t(nall);
257 //
258 if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1;
259
260
261}
262
263//
264//
265void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed, AliESD */*event*/)
266{
267 //
268 //
269 //calculates derived variables
270 //
271 //
272 UpdatePoints(fESDtrack);
273 fBestTOFmatch=1000;
274 AliTrackReference * ref = &(info->fTrackRef);
275 fTPCinR0[0] = info->fTrackRef.X();
276 fTPCinR0[1] = info->fTrackRef.Y();
277 fTPCinR0[2] = info->fTrackRef.Z();
278 fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
279 fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
280 //
281 fTPCinP0[0] = ref->Px();
282 fTPCinP0[1] = ref->Py();
283 fTPCinP0[2] = ref->Pz();
284 fTPCinP0[3] = ref->Pt();
285 fTPCinP0[4] = ref->P();
286 fDeltaP = (ref->P()-info->fParticle.P())/info->fParticle.P();
287 //
288 //
289 if (fTPCinP0[3]>0.0000001){
290 //
291 fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
292 fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
293 }
294 //
295 //
296 fITSinP0[0]=info->fParticle.Px();
297 fITSinP0[1]=info->fParticle.Py();
298 fITSinP0[2]=info->fParticle.Pz();
299 fITSinP0[3]=info->fParticle.Pt();
300 //
301 fITSinR0[0]=info->fParticle.Vx();
302 fITSinR0[1]=info->fParticle.Vy();
303 fITSinR0[2]=info->fParticle.Vz();
304 fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
305 fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
306 //
307 //
308 if (fITSinP0[3]>0.0000001){
309 fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
310 fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
311 }
312 //
313 for (Int_t i=0;i<4;i++) fStatus[i] =0;
314 fReconstructed = kFALSE;
315 fTPCOn = kFALSE;
316 fITSOn = kFALSE;
317 fTRDOn = kFALSE;
318 if (reconstructed==kFALSE) return;
319
320 fLabels[0] = info->fLabel;
321 fLabels[1] = info->fPrimPart;
322 fReconstructed = kTRUE;
323 fTPCOn = ((fESDtrack->GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
324 fITSOn = ((fESDtrack->GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
325 fTRDOn = ((fESDtrack->GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
326 //
327 //
328 if ((fESDtrack->GetStatus()&AliESDtrack::kTPCrefit)>0){
329 fStatus[1] =3;
330 }
331 else{
332 if ((fESDtrack->GetStatus()&AliESDtrack::kTPCout)>0){
333 fStatus[1] =2;
334 }
335 else{
336 if ((fESDtrack->GetStatus()&AliESDtrack::kTPCin)>0)
337 fStatus[1]=1;
338 }
339 }
340 //
341 if ((fESDtrack->GetStatus()&AliESDtrack::kITSout)>0){
342 fStatus[0] =2;
343 }
344 else{
345 if ((fESDtrack->GetStatus()&AliESDtrack::kITSrefit)>0){
346 fStatus[0] =1;
347 }
348 else{
349 fStatus[0]=0;
350 }
351 }
352
353 //
354 //
355 if ((fESDtrack->GetStatus()&AliESDtrack::kTRDrefit)>0){
356 fStatus[2] =2;
357 }
358 else{
359 if ((fESDtrack->GetStatus()&AliESDtrack::kTRDout)>0){
360 fStatus[2] =1;
361 }
362 }
363 if ((fESDtrack->GetStatus()&AliESDtrack::kTRDStop)>0){
364 fStatus[2] =10;
365 }
366
367 //
368 //TOF
369 //
370 if (((fESDtrack->GetStatus()&AliESDtrack::kTOFout)>0)){
371 //
372 // best tof match
373 Double_t times[5];
374 fESDtrack->GetIntegratedTimes(times);
375 for (Int_t i=0;i<5;i++){
376 if ( TMath::Abs(fESDtrack->GetTOFsignal()-times[i]) <TMath::Abs(fBestTOFmatch) ){
377 fBestTOFmatch = fESDtrack->GetTOFsignal()-times[i];
378 }
379 }
380 Int_t toflabel[3];
381 fESDtrack->GetTOFLabel(toflabel);
382 Bool_t toffake=kTRUE;
383 Bool_t tofdaughter=kFALSE;
384 for (Int_t i=0;i<3;i++){
385 if (toflabel[i]<0) continue;
386 if (toflabel[i]== TMath::Abs(fESDtrack->GetLabel())) toffake=kFALSE;
387 if (toflabel[i]==info->fParticle.GetDaughter(0) || (toflabel[i]==info->fParticle.GetDaughter(1))) tofdaughter=kTRUE; // decay product of original particle
388 fStatus[3]=1;
389 }
390 if (toffake) fStatus[3] =3; //total fake
391 if (tofdaughter) fStatus[3]=2; //fake because of decay
392 }else{
393 fStatus[3]=0;
394 }
395
396
397 if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){
398 //TPC
399 fESDtrack->GetInnerXYZ(fTPCinR1);
400 fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
401 fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);
402 fESDtrack->GetInnerPxPyPz(fTPCinP1);
403 fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
404 fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]);
405 //
406 //
407 if (fTPCinP1[3]>0.000000000000001){
408 fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
409 fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);
410 }
411 Double_t cov[15], param[5],x, alpha;
412 fESDtrack->GetInnerExternalCovariance(cov);
413 fESDtrack->GetInnerExternalParameters(alpha, x,param);
414 if (x<50) return ;
415 //
416 fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3]; //delta rfi
417 fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
418 fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]); //delta z
419 fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
420 fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
421 fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
422 fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]));
423 fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
424 fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
425 Double_t sign = (param[4]>0)? 1.:-1;
426 fSign =sign;
427 fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(TMath::Abs(cov[14]));
428 }
429 if (fITSOn){
430 // ITS
431 Double_t param[5],x;
432 fESDtrack->GetExternalParameters(x,param);
433 // fESDtrack->GetConstrainedExternalParameters(x,param);
434 Double_t cov[15];
435 fESDtrack->GetExternalCovariance(cov);
436 //fESDtrack->GetConstrainedExternalCovariance(cov);
437 if (TMath::Abs(param[4])<0.0000000001) return;
438
439 fESDtrack->GetXYZ(fITSinR1);
440 fESDtrack->GetPxPyPz(fITSinP1);
441 fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]);
442 //
443 fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]);
444 fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]);
445 //
446 //
447 if (fITSinP1[3]>0.0000001){
448 fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]);
449 fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]);
450 }
451 //
452 //
453 fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3]; //delta rfi
454 fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]);
455 fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]); //delta z
456 fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]);
457 fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]);
458 fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]);
459 fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1]));
460 fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
461 fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);
462 Double_t sign = (param[4]>0) ? 1:-1;
463 fSign = sign;
464 fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);
465 }
466
467}
468
469
470void AliESDRecV0Info::Update(Float_t vertex[3])
471{
472
473 if ( (fT1.fStatus[1]>0)&& (fT2.fStatus[1]>0)){
474 Float_t distance1,distance2;
475 Double_t xx[3],pp[3];
476 //
477 Double_t xd[3],pd[3],signd;
478 Double_t xm[3],pm[3],signm;
479 //
480 //
481 if (fT1.fITSOn&&fT2.fITSOn){
482 for (Int_t i=0;i<3;i++){
483 xd[i] = fT2.fITSinR1[i];
484 pd[i] = fT2.fITSinP1[i];
485 xm[i] = fT1.fITSinR1[i];
486 pm[i] = fT1.fITSinP1[i];
487 }
488 }
489 else{
490
491 for (Int_t i=0;i<3;i++){
492 xd[i] = fT2.fTPCinR1[i];
493 pd[i] = fT2.fTPCinP1[i];
494 xm[i] = fT1.fTPCinR1[i];
495 pm[i] = fT1.fTPCinP1[i];
496 }
497 }
498 //
499 //
500 signd = fT2.fSign<0 ? -1:1;
501 signm = fT1.fSign<0 ? -1:1;
502
503 AliHelix dhelix1(xd,pd,signd);
504 dhelix1.GetMomentum(0,pp,0);
505 dhelix1.Evaluate(0,xx);
506 //
507 // Double_t x2[3],p2[3];
508 //
509 AliHelix mhelix(xm,pm,signm);
510 //
511 //find intersection linear
512 //
513 Double_t phase[2][2],radius[2];
514 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
515 Double_t delta1=10000,delta2=10000;
516
517 if (points==1){
518 fRs[0] = TMath::Sqrt(radius[0]);
519 fRs[1] = TMath::Sqrt(radius[0]);
520 }
521 if (points==2){
522 fRs[0] =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
523 fRs[1] =TMath::Max(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
524 }
525
526 if (points>0){
527 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
528 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
529 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
530 }
531 if (points==2){
532 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
533 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
534 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
535 }
536 if (points==1){
537 fRs[0] = TMath::Sqrt(radius[0]);
538 fRs[1] = TMath::Sqrt(radius[0]);
539 fDistMinR = delta1;
540 }
541 if (points==2){
542 if (radius[0]<radius[1]){
543 fRs[0] = TMath::Sqrt(radius[0]);
544 fRs[1] = TMath::Sqrt(radius[1]);
545 fDistMinR = delta1;
546 }
547 else{
548 fRs[0] = TMath::Sqrt(radius[1]);
549 fRs[1] = TMath::Sqrt(radius[0]);
550 fDistMinR = delta2;
551 }
552 }
553 //
554 //
555 distance1 = TMath::Min(delta1,delta2);
556 //
557 //find intersection parabolic
558 //
559 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
560 delta1=10000,delta2=10000;
561
562 if (points>0){
563 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
564 }
565 if (points==2){
566 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
567 }
568
569 distance2 = TMath::Min(delta1,delta2);
570 if (distance2>100) fDist2 =100;
571 return;
572 if (delta1<delta2){
573 //get V0 info
574 dhelix1.Evaluate(phase[0][0],fXr);
575 dhelix1.GetMomentum(phase[0][0],fPdr);
576 mhelix.GetMomentum(phase[0][1],fPm);
577 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
578 fRr = TMath::Sqrt(radius[0]);
579 }
580 else{
581 dhelix1.Evaluate(phase[1][0],fXr);
582 dhelix1.GetMomentum(phase[1][0], fPdr);
583 mhelix.GetMomentum(phase[1][1], fPm);
584 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
585 fRr = TMath::Sqrt(radius[1]);
586 }
587 fDist1 = TMath::Sqrt(distance1);
588 fDist2 = TMath::Sqrt(distance2);
589
590 if (fDist2<10.5){
591 Double_t x,alpha,param[5],cov[15];
592 //
593 fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
594 fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
595 AliExternalTrackParam paramm(x,alpha,param,cov);
596 //
597 fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
598 fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
599 AliExternalTrackParam paramd(x,alpha,param,cov);
600 }
601 //
602 //
603
604 Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
605 Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
606
607 Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
608 Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
609 vnorm2 = TMath::Sqrt(vnorm2);
610 Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
611 Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
612 pnorm2 = TMath::Sqrt(pnorm2);
613
614 fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
615 fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
616 fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
617 }
618}
619
620////
621void AliESDRecKinkInfo::Update()
622{
623
624 if ( (fT1.fTPCOn)&& (fT2.fTPCOn)){
625 //
626 // IF BOTH RECONSTRUCTED
627 Float_t distance1,distance2;
628 Double_t xx[3],pp[3];
629 //
630 Double_t xd[3],pd[3],signd;
631 Double_t xm[3],pm[3],signm;
632 for (Int_t i=0;i<3;i++){
633 xd[i] = fT2.fTPCinR1[i];
634 pd[i] = fT2.fTPCinP1[i];
635 xm[i] = fT1.fTPCinR1[i];
636 pm[i] = fT1.fTPCinP1[i];
637 }
638 signd = fT2.fSign<0 ? -1:1;
639 signm = fT1.fSign<0 ? -1:1;
640
641 AliHelix dhelix1(xd,pd,signd);
642 dhelix1.GetMomentum(0,pp,0);
643 dhelix1.Evaluate(0,xx);
644 //
645 // Double_t x2[3],p2[3];
646 //
647 AliHelix mhelix(xm,pm,signm);
648 //
649 //find intersection linear
650 //
651 Double_t phase[2][2],radius[2];
652 Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
653 Double_t delta1=10000,delta2=10000;
654
655 if (points==1){
656 fMinR = TMath::Sqrt(radius[0]);
657 }
658 if (points==2){
659 fMinR =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
660 }
661
662 if (points>0){
663 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
664 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
665 dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
666 }
667 if (points==2){
668 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
669 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
670 dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
671 }
672 if (points==1){
673 fMinR = TMath::Sqrt(radius[0]);
674 fDistMinR = delta1;
675 }
676 if (points==2){
677 if (radius[0]<radius[1]){
678 fMinR = TMath::Sqrt(radius[0]);
679 fDistMinR = delta1;
680 }
681 else{
682 fMinR = TMath::Sqrt(radius[1]);
683 fDistMinR = delta2;
684 }
685 }
686 //
687 //
688 distance1 = TMath::Min(delta1,delta2);
689 //
690 //find intersection parabolic
691 //
692 points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
693 delta1=10000,delta2=10000;
694
695 if (points>0){
696 dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
697 }
698 if (points==2){
699 dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
700 }
701
702 distance2 = TMath::Min(delta1,delta2);
703 if (delta1<delta2){
704 //get V0 info
705 dhelix1.Evaluate(phase[0][0],fXr);
706 dhelix1.GetMomentum(phase[0][0],fPdr);
707 mhelix.GetMomentum(phase[0][1],fPm);
708 dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
709 fRr = TMath::Sqrt(radius[0]);
710 }
711 else{
712 dhelix1.Evaluate(phase[1][0],fXr);
713 dhelix1.GetMomentum(phase[1][0], fPdr);
714 mhelix.GetMomentum(phase[1][1], fPm);
715 dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
716 fRr = TMath::Sqrt(radius[1]);
717 }
718 fDist1 = TMath::Sqrt(distance1);
719 fDist2 = TMath::Sqrt(distance2);
720
721 if (fDist2<10.5){
722 Double_t x,alpha,param[5],cov[15];
723 //
724 fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
725 fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
726 AliExternalTrackParam paramm(x,alpha,param,cov);
727 //
728 fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
729 fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
730 AliExternalTrackParam paramd(x,alpha,param,cov);
731 /*
732 AliESDkink kink;
733 kink.Update(&paramm,&paramd);
734 // kink.Dump();
735 Double_t diff = kink.fRr-fRr;
736 Double_t diff2 = kink.fDist2-fDist2;
737 printf("Diff\t%f\t%f\n",diff,diff2);
738 */
739 }
740
741 //
742 //
743 }
744
745}
746