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 #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
  87 ClassImp(AliESDRecInfo)
  88 ClassImp(AliESDRecV0Info)
  89 ClassImp(AliESDRecKinkInfo)
  90
  91
  92
  93
  94 AliTPCParam * GetTPCParam(){
  95   AliTPCParamSR * par = new AliTPCParamSR;
  96   par->Update();
  97   return par;
  98 }
  99
 100
 101
 102
 103 AliESDRecInfo::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
 125 AliESDRecInfo::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
 137 AliESDRecInfo::~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
 153 void 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
 168 void 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
 185 void  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 //
 265 void 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
 470 void  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 ////
 621 void  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