PWG1/AliRecInfo.cxx

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