Adding new code for kink and V0 reconstruction (AliESDkink and AliESDV0MI classes...
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 15 Oct 2004 14:48:52 +0000 (14:48 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 15 Oct 2004 14:48:52 +0000 (14:48 +0000)
30 files changed:
ITS/AliITStrackerMI.cxx
ITS/AliITStrackerMI.h
ITS/ITSLinkDef.h
STEER/AliESD.cxx
STEER/AliESD.h
STEER/AliESDComparisonMI.C
STEER/AliESDComparisonMI.h
STEER/AliESDV0MI.cxx [new file with mode: 0644]
STEER/AliESDV0MI.h [new file with mode: 0644]
STEER/AliESDkink.cxx [new file with mode: 0644]
STEER/AliESDkink.h [new file with mode: 0644]
STEER/AliESDtrack.cxx
STEER/AliESDtrack.h
STEER/AliExternalTrackParam.cxx [new file with mode: 0644]
STEER/AliExternalTrackParam.h [new file with mode: 0644]
STEER/AliGenInfo.C
STEER/AliGenInfo.h
STEER/AliHelix.cxx [moved from TPC/AliHelix.cxx with 90% similarity]
STEER/AliHelix.h [moved from TPC/AliHelix.h with 97% similarity]
STEER/AliTrackParam.cxx [new file with mode: 0644]
STEER/AliTrackParam.h [new file with mode: 0644]
STEER/ESDLinkDef.h
STEER/STEERLinkDef.h
STEER/libESD.pkg
TPC/AliTPCtrack.cxx
TPC/AliTPCtrack.h
TPC/AliTPCtrackerMI.cxx
TPC/AliTPCtrackerMI.h
TPC/TPCrecLinkDef.h
TPC/libTPCrec.pkg

index c44c961aacebad29ee4997777a3f758b57b3bf80..47af8e6863bcb842bfc626a3758f2b79f5d47067 100644 (file)
 #include "AliITSclusterV2.h"
 #include "AliITStrackerMI.h"
 #include "TMatrixD.h"
+#include "TFile.h"
+#include "TTree.h"
 #include "AliHelix.h"
-
+#include "AliESDV0MI.h"
 
 
 ClassImp(AliITStrackerMI)
-ClassImp(AliITSRecV0Info)
 
 
 
@@ -227,7 +228,7 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
       if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
       if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
       if (esd->GetStatus()&AliESDtrack::kITSin) continue;
-
+      if (esd->GetKinkIndex(0)>0) continue;   //kink daughter
       AliITStrackMI *t=0;
       try {
         t=new AliITStrackMI(*esd);
@@ -319,7 +320,7 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
   }
 
   //GetBestHypothesysMIP(itsTracks);
-  //FindV0(event);
+  FindV0(event);
 
   itsTracks.Delete();
   //
@@ -3012,15 +3013,25 @@ void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
 
 
 
-void  AliITStrackerMI::FindV0(AliESD */*event*/)
+
+
+
+
+void  AliITStrackerMI::FindV0(AliESD *event)
 {
   //
   // fast V0 finder
   //
+  //fV0Array->Clean();
   AliHelix helixes[30000];
   TObjArray trackarray(30000);
   Float_t dist[30000];
-  AliITSRecV0Info vertex;
+  Float_t mindist[30000];
+  AliESDV0MI  *vertexarray    = new AliESDV0MI[20000];
+  AliV0vertex *oldvertexarray = new AliV0vertex[20000];
+  AliESDV0MI *pvertex     = &vertexarray[0];
+  
+
   //
   Int_t entries = fTrackHypothesys.GetEntriesFast();
   for (Int_t i=0;i<entries;i++){
@@ -3029,44 +3040,171 @@ void  AliITStrackerMI::FindV0(AliESD */*event*/)
     AliITStrackMI * track = (AliITStrackMI*)array->At(fBestTrackIndex[i]);
     if (track){
       dist[i] = TMath::Sqrt(track->fD[0]*track->fD[0]+track->fD[1]*track->fD[1]);
+      Float_t pt1 = TMath::Abs(track->fP4*track->GetConvConst());
+      mindist[i] = TMath::Min(0.02+0.03*pt1,0.15);
+      if (mindist[i]<0.05) mindist[i]=0.05;
       trackarray.AddAt(track,i);
       new (&helixes[i]) AliHelix(*track);
     }
   }
+  //  Int_t multifound=0;
+  Int_t vertexall =0;
+  AliESDV0MI tempvertex;
+  Float_t primvertex[3]={GetX(),GetY(),GetZ()};
+      
   for (Int_t itrack0=0;itrack0<entries;itrack0++){
     AliITStrackMI * track0 = (AliITStrackMI*)trackarray.At(itrack0);
     if (!track0) continue;
-    if (dist[itrack0]<0.2) continue;
-    for (Int_t itrack1=itrack0+1;itrack1<entries;itrack1++){
+    if (track0->fP4>0) continue;
+    if (dist[itrack0]<mindist[itrack0]) continue;
+    Int_t vertexes =0;
+    //
+    for (Int_t itrack1=0;itrack1<entries;itrack1++){
       AliITStrackMI * track1 = (AliITStrackMI*)trackarray.At(itrack1);
       if (!track1) continue;
-      if (dist[itrack1]<0.2) continue;
-      if (track1->fP4*track0->fP4>0) continue; //the same sign
+      if (track1->fP4<0) continue;
+      if (dist[itrack1]<mindist[itrack1]) continue;
+      //
       AliHelix *h1 = &helixes[itrack0];
-      AliHelix *h2 = &helixes[itrack1];
-      
-      Double_t distance = TestV0(h1,h2,&vertex);
+      AliHelix *h2 = &helixes[itrack1];      
+      Double_t distance = TestV0(h1,h2,pvertex);
       if (distance>0.4) continue;
-      if (vertex.fRr<0.3) continue;
-      if (vertex.fRr>27) continue;
-      Float_t v[3]={GetX(),GetY(),GetZ()};
-      vertex.Update(v,track0->fMass, track1->fMass);
+      //
+      if (distance>0.3*(dist[itrack1]+dist[itrack0])) continue;
+      //if (distance>0.2*dist[itrack0]) continue;
+      if (pvertex->fRr<0.3*(dist[itrack1]+dist[itrack0])) continue;
+      if (pvertex->fRr>27) continue;
+      pvertex->SetM(*track0);
+      pvertex->SetP(*track1);
+      pvertex->Update(primvertex);
+      //
+      //
+      if (pvertex->fPointAngle<0.85) continue;       
+      //
+      //
+      pvertex->fLab[0]=track0->GetLabel();
+      pvertex->fLab[1]=track1->GetLabel();
+      pvertex->fIndex[0] = track0->GetESDtrack()->GetID();
+      pvertex->fIndex[1] = track1->GetESDtrack()->GetID();
+      // calculate chi2s
+      //
+      pvertex->fChi2After  = 0;
+      pvertex->fChi2Before = 0;
+      pvertex->fNBefore=0;
+      pvertex->fNAfter=0;
       
-      if ((TMath::Abs(vertex.fInvMass-0.4976)<0.05 || TMath::Abs(vertex.fInvMass-0.0)<0.1||  TMath::Abs(vertex.fInvMass-1.1)<0.1)) {
-       if (vertex.fPointAngle<0.85) continue;  
+      for (Int_t i=0;i<6;i++){
+       Float_t radius = fgLayers[i].GetR();
+       if (pvertex->fRr>radius+0.5){
+         pvertex->fNBefore+=2.;
+         //
+         if (track0->fClIndex[i]<=0) {
+           pvertex->fChi2Before+=9;
+         }else{
+           Float_t chi2 = track0->fDy[i]*track0->fDy[i]/(track0->fSigmaY[i]*track0->fSigmaY[i])+
+             track0->fDz[i]*track0->fDz[i]/(track0->fSigmaZ[i]*track0->fSigmaZ[i]);
+             pvertex->fChi2Before+=chi2;
+         }
+
+         if (track1->fClIndex[i]<=0) {
+           pvertex->fChi2Before+=9;
+         }else{
+           Float_t chi2 = track1->fDy[i]*track1->fDy[i]/(track1->fSigmaY[i]*track1->fSigmaY[i])+
+             track1->fDz[i]*track1->fDz[i]/(track1->fSigmaZ[i]*track1->fSigmaZ[i]);
+             pvertex->fChi2Before+=chi2;
+         }
+       }
+
+       if (pvertex->fRr<radius-0.5){
+         pvertex->fNAfter+=2.;
+         //
+         if (track0->fClIndex[i]<=0) {
+           pvertex->fChi2After+=9;
+         }else{
+           Float_t chi2 = track0->fDy[i]*track0->fDy[i]/(track0->fSigmaY[i]*track0->fSigmaY[i])+
+             track0->fDz[i]*track0->fDz[i]/(track0->fSigmaZ[i]*track0->fSigmaZ[i]);
+             pvertex->fChi2After+=chi2;
+         }
+
+         if (track1->fClIndex[i]<=0) {
+           pvertex->fChi2After+=9;
+         }else{
+           Float_t chi2 = track1->fDy[i]*track1->fDy[i]/(track1->fSigmaY[i]*track1->fSigmaY[i])+
+             track1->fDz[i]*track1->fDz[i]/(track1->fSigmaZ[i]*track1->fSigmaZ[i]);
+             pvertex->fChi2After+=chi2;
+         }
+       }
       }
-      else{
-       if (vertex.fPointAngle<0.98) continue;  
+      if (pvertex->fNBefore>2){
+       if (pvertex->fChi2Before/pvertex->fNBefore<4.) continue; //have clusters before vetex
+      }
+      distance = FindBestPair(itrack0,itrack1,pvertex);
+      if (pvertex->fPointAngle<0.85) continue; 
+      //
+      if (distance>0.3) continue;
+      //      if (pvertex->fDistSigma*6>pvertex->fRr) continue;
+      // if (pvertex->fDistSigma>0.4) continue;
+      //if (pvertex->fDistNorm>5.5) continue;
+      new (&oldvertexarray[vertexall]) AliV0vertex(*track0,*track1) ;
+      vertexes++;
+      vertexall++;
+      pvertex = &vertexarray[vertexall];
+    }
+  }
+  //  printf("\n\n\nMultifound\t%d\n\n\n",multifound);
+  //
+  // sort vertices according quality
+  Float_t quality[10000];
+  Int_t   indexes[10000];
+  Int_t   trackvertices[30000];
+  for (Int_t i=0;i<entries;i++) trackvertices[i]=0;
+  for (Int_t i=0;i<vertexall;i++) quality[i] = 1.-vertexarray[i].fPointAngle;
+  TMath::Sort(vertexall,quality,indexes,kFALSE);
+  
+  for (Int_t i=0;i<vertexall;i++){
+    pvertex= &vertexarray[indexes[i]];
+    Int_t index0 = vertexarray[indexes[i]].fIndex[0];
+    Int_t index1 = vertexarray[indexes[i]].fIndex[1];
+    vertexarray[indexes[i]].fOrder[2] = i;
+    vertexarray[indexes[i]].fOrder[1] = trackvertices[index1];
+    vertexarray[indexes[i]].fOrder[0] = trackvertices[index0];
+    if (trackvertices[index1]+trackvertices[index0]>2) continue;
+    if (trackvertices[index1]+trackvertices[index0]>0) {
+      //      if (pvertex->fPointAngle<0.995)  continue;
+    }
+    trackvertices[index0]++;
+    trackvertices[index1]++;
+    
+    //    event->AddV0(&oldvertexarray[indexes[i]]);
+    Int_t v0index = event->AddV0MI(&vertexarray[indexes[i]]);
+    AliESDtrack * ptrack0 = event->GetTrack(vertexarray[indexes[i]].fIndex[0]);
+    AliESDtrack * ptrack1 = event->GetTrack(vertexarray[indexes[i]].fIndex[1]);
+    if (!ptrack0 || !ptrack1){
+      printf("BBBBBBBUUUUUUUUUUGGGGGGGGGG\n");
+    }
+    Int_t v0index0[3]={ptrack0->GetV0Index(0),ptrack0->GetV0Index(1),ptrack0->GetV0Index(2)};
+    Int_t v0index1[3]={ptrack1->GetV0Index(0),ptrack1->GetV0Index(1),ptrack1->GetV0Index(2)};
+    for (Int_t i=0;i<3;i++){
+      if (v0index0[i]<0) {
+       v0index0[i]=v0index;
+       ptrack0->SetV0Indexes(v0index0);
+       break;
       }
-      if (TMath::Abs(TMath::Abs(track0->fLab)-TMath::Abs(track1->fLab))<2){
-       Float_t mindist = FindBestPair(itrack0,itrack1,&vertex);
-       if (mindist>1) FindBestPair(itrack0,itrack1,&vertex);
+    }
+    for (Int_t i=0;i<3;i++){
+      if (v0index1[i]<0) {
+       v0index1[i]=v0index;
+       ptrack1->SetV0Indexes(v0index1);
+       break;
       }
-    }    
+    }
   }
+  delete[] vertexarray;
+  delete[] oldvertexarray;
 }
 
-Double_t  AliITStrackerMI::FindBestPair(Int_t esdtrack0, Int_t esdtrack1,AliITSRecV0Info *vertex)
+
+Double_t  AliITStrackerMI::FindBestPair(Int_t esdtrack0, Int_t esdtrack1, AliESDV0MI *vertex)
 {
   //
   // try to find best pair from the tree of track hyp.
@@ -3075,89 +3213,67 @@ Double_t  AliITStrackerMI::FindBestPair(Int_t esdtrack0, Int_t esdtrack1,AliITSR
   Int_t entries0 = array0->GetEntriesFast();
   TObjArray * array1 = (TObjArray*)fTrackHypothesys.At(esdtrack1);
   Int_t entries1 = array1->GetEntriesFast();  
+  //  Float_t primvertex[3]={GetX(),GetY(),GetZ()};
   //
   //
-  AliITSRecV0Info v0;
-  AliITSRecV0Info vbest;
+  //AliESDV0MI v0;
   Double_t criticalradius = vertex->fRr;
-  Double_t mindist =1000;
-  Int_t bestpair[2];
   //
+  AliITStrackMI * track0= (AliITStrackMI*)array0->At(fBestTrackIndex[esdtrack0]);
+  AliITStrackMI * track1= (AliITStrackMI*)array1->At(fBestTrackIndex[esdtrack1]);
+  //
+  // find the best tracks after decay point
   for (Int_t itrack0=0;itrack0<entries0;itrack0++){
-    AliITStrackMI * track0 = (AliITStrackMI*)array0->At(itrack0);
-    if (!track0) continue;
-    if (track0->fX<criticalradius-1) continue;
-    if (track0->fX>criticalradius+5) continue;
-    for (Int_t itrack1=0;itrack1<entries1;itrack1++){
-      AliITStrackMI * track1 = (AliITStrackMI*)array1->At(itrack1);
-      if (!track1) continue;      
-      if (track1->fX<criticalradius-1) continue;
-      if (track1->fX>criticalradius+5) continue;
-
-      AliHelix h0(*track0);
-      AliHelix h1(*track1);
-      Double_t distance = TestV0(&h0,&h1,&v0);
-      if (v0.fRr>30) continue;
-      if (v0.fRr<0.2) continue;
-      Float_t v[3]={GetX(),GetY(),GetZ()};
-      v0.Update(v,track0->fMass, track1->fMass);
-      if (distance<mindist){
-       mindist = distance;
-       bestpair[0]=itrack0;
-       bestpair[1]=itrack1;
-       vbest = v0;     
-      }      
+    AliITStrackMI * track = (AliITStrackMI*)array0->At(itrack0);
+    if (!track) continue;
+    if (track->fX<criticalradius-1) continue;
+    if (track->fX>criticalradius){
+      track0 = track;
+      break;  
     }
   }
-  if (mindist<0.2){
-    vbest.Dump();
-    vertex->Dump();
-  }
-  return mindist;
-}
-
-
-void  AliITSRecV0Info::Update(Float_t vertex[3], Float_t mass1, Float_t mass2)
-{
-  //
-  //Update V0 information
-  //
-  Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
-  Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
-
-  Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
-  Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
-  vnorm2 = TMath::Sqrt(vnorm2);
-  Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
-  Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
-  pnorm2 = TMath::Sqrt(pnorm2);
-  
-  fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
-  fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
-  fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
-  
-  Float_t e1    = TMath::Sqrt(mass1*mass1+
-                             fPdr[0]*fPdr[0]+
-                             fPdr[1]*fPdr[1]+
-                             fPdr[2]*fPdr[2]);
-  Float_t e2    = TMath::Sqrt(mass2*mass2+
-                             fPm[0]*fPm[0]+
-                             fPm[1]*fPm[1]+
-                             fPm[2]*fPm[2]);
-  
-  fInvMass =  
-    (fPm[0]+fPdr[0])*(fPm[0]+fPdr[0])+
-    (fPm[1]+fPdr[1])*(fPm[1]+fPdr[1])+
-    (fPm[2]+fPdr[2])*(fPm[2]+fPdr[2]);
-  
-  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
-
 
+  for (Int_t itrack1=0;itrack1<entries1;itrack1++){
+    AliITStrackMI * track = (AliITStrackMI*)array1->At(itrack1);
+    if (!track) continue;
+    if (track->fX<criticalradius-1) continue;
+    if (track->fX>criticalradius){
+      track1 = track;
+      break;  
+    }
+  }
+  //propagate to vertex
+  Double_t alpha = TMath::ATan2(vertex->fXr[1],vertex->fXr[0]);
+  Double_t radius =TMath::Sqrt(vertex->fXr[1]*vertex->fXr[1]+vertex->fXr[0]*vertex->fXr[0]);
+  AliITStrackMI track0p = *track0;
+  AliITStrackMI track1p = *track1;
+  if (!track0p.Propagate(alpha,radius)) return 100000;
+  if (!track1p.Propagate(alpha,radius)) return 100000;
+  //track0p.Propagate(alpha,radius);
+  //track1p.Propagate(alpha,radius);
+  //
+  //
+  //vertex2->SetM(track0p);
+  //vertex2->SetP(track1p);
+  //vertex2->Update(primvertex);
+  //
+  AliHelix h0(track0p);
+  AliHelix h1(track1p);
+  Double_t distance = TestV0(&h0,&h1,vertex);
+  Float_t v[3]={GetX(),GetY(),GetZ()};
+  //  vertex->Update(v,track0->fMass, track1->fMass);  
+  vertex->SetM(*track0);
+  vertex->SetP(*track1);
+  vertex->Update(v);
+  Float_t sigma = TMath::Sqrt(track1p.GetSigmaY2()+track1p.GetSigmaZ2()+track0p.GetSigmaY2()+track0p.GetSigmaZ2());
+  vertex->fDistNorm = distance/sigma;
+  vertex->fDistSigma = sigma;
+  return distance;
 }
 
 
 
-Double_t   AliITStrackerMI::TestV0(AliHelix *helix1, AliHelix *helix2, AliITSRecV0Info *vertex)
+Double_t   AliITStrackerMI::TestV0(AliHelix *helix1, AliHelix *helix2, AliESDV0MI *vertex)
 {
   //
   // test the helixes for the distnce calculate vertex characteristic
@@ -3206,15 +3322,15 @@ Double_t   AliITStrackerMI::TestV0(AliHelix *helix1, AliHelix *helix2, AliITSRec
     if (delta1<delta2){
       //get V0 info
       dhelix1.Evaluate(phase[0][0],vertex->fXr);
-      dhelix1.GetMomentum(phase[0][0],vertex->fPdr);
-      mhelix.GetMomentum(phase[0][1],vertex->fPm);
+      dhelix1.GetMomentum(phase[0][0],vertex->fPP);
+      mhelix.GetMomentum(phase[0][1],vertex->fPM);
       dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],vertex->fAngle);
        vertex->fRr = TMath::Sqrt(radius[0]);
     }
     else{
       dhelix1.Evaluate(phase[1][0],vertex->fXr);
-      dhelix1.GetMomentum(phase[1][0],vertex->fPdr);
-      mhelix.GetMomentum(phase[1][1],vertex->fPm);
+      dhelix1.GetMomentum(phase[1][0],vertex->fPP);
+      mhelix.GetMomentum(phase[1][1],vertex->fPM);
       dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],vertex->fAngle);
       vertex->fRr = TMath::Sqrt(radius[1]);
     }
@@ -3224,3 +3340,9 @@ Double_t   AliITStrackerMI::TestV0(AliHelix *helix1, AliHelix *helix2, AliITSRec
   return  vertex->fDist2;
 }
 
+
+
+
+
+
+
index c58f7b1dc7c68bbdf23bb97e60332ae1af2dc042..3dcd7bd21ead829ae44f300b7c019a800a708f6f 100644 (file)
 #include "AliTracker.h"
 #include "AliITStrackMI.h"
 #include "AliITSclusterV2.h"
+#include "AliV0vertex.h"
 
 class AliESD;
 class AliITSgeom;
 class TTree;
 class AliHelix;
 class AliV0vertex;
-
-
-class AliITSRecV0Info: public TObject {
-  friend class AliITStrackerMI;
-protected:
-  void Update(Float_t vertex[3], Float_t mass1, Float_t mass2);
-  Double_t       fDist1;    //info about closest distance according closest MC - linear DCA
-  Double_t       fDist2;    //info about closest distance parabolic DCA
-  Double_t       fInvMass;  //reconstructed invariant mass -
-  //
-  Double_t       fPdr[3];    //momentum at vertex daughter  - according approx at DCA
-  Double_t       fXr[3];     //rec. position according helix
-  //
-  Double_t       fPm[3];    //momentum at the vertex mother
-  Double_t       fAngle[3]; //three angles
-  Double_t       fRr;       // rec position of the vertex 
-  Int_t          fLab[2];   //MC label of the partecle 
-  Float_t        fPointAngleFi; //point angle fi
-  Float_t        fPointAngleTh; //point angle theta
-  Float_t        fPointAngle;   //point angle full
-  ClassDef(AliITSRecV0Info,1)  // container for  
-};
+class AliESDV0MI;
 
 
 
@@ -190,8 +170,8 @@ public:
 
 protected:
   void FindV0(AliESD *event);  //try to find V0
-  Double_t  TestV0(AliHelix *h1, AliHelix *h2, AliITSRecV0Info *vertex);  //try to find V0 - return DCA
-  Double_t  FindBestPair(Int_t esdtrack0, Int_t esdtrack1,AliITSRecV0Info *vertex);  // try to find best pair from the tree of track hyp.
+  Double_t  TestV0(AliHelix *h1, AliHelix *h2, AliESDV0MI *vertex);  //try to find V0 - return DCA
+  Double_t  FindBestPair(Int_t esdtrack0, Int_t esdtrack1,AliESDV0MI *vertex);  // try to find best pair from the tree of track hyp.
   void CookLabel(AliKalmanTrack *t,Float_t wrong) const;
   void CookLabel(AliITStrackMI *t,Float_t wrong) const;
   Double_t GetEffectiveThickness(Double_t y, Double_t z) const;
@@ -235,6 +215,7 @@ protected:
   AliITStrackMI fBestTrack;              // "best" track 
   AliITStrackMI fTrackToFollow;          // followed track
   TObjArray     fTrackHypothesys;        // ! array with track hypothesys- ARRAY is the owner of tracks- MI
+  TObjArray     fV0Array;                // ! array of V0
   Int_t         fBestTrackIndex[100000]; // ! index of the best track
   Int_t         fCurrentEsdTrack;        // ! current esd track           - MI
   Int_t fPass;                           // current pass through the data 
index 8058b89d66479a213993f20873a1c6c56c461934..dffd0ba209d96851a53065840b051d0f03dd3bf1 100644 (file)
 #pragma link C++ class AliITStrackerV2+;
 #pragma link C++ class AliITStrackMI+;
 #pragma link C++ class AliITStrackerMI+;
-#pragma link C++ class AliITSRecV0Info+;
+//#pragma link C++ class AliITSRecV0Info+;
 #pragma link C++ class  AliV0vertex+;
 #pragma link C++ class  AliV0vertexer+;
 #pragma link C++ class  AliCascadeVertex+;
index 5f1aaaa390fc34d1844ddbddebd1d3f8580b0c55..55230c579cda5c820c15d00b2873abe6e7d0c644 100644 (file)
@@ -44,8 +44,10 @@ AliESD::AliESD():
   fHLTHoughTracks("AliESDHLTtrack",15000),
   fMuonTracks("AliESDMuonTrack",30),
   fPmdTracks("AliESDPmdTrack",3000),
-  fV0s("AliESDv0",200),
+  fV0s("AliESDv0",200),  
   fCascades("AliESDcascade",20),
+  fKinks("AliESDkink",4000),
+  fV0MIs("AliESDV0MI",4000),
   fPHOSParticles(0), 
   fEMCALParticles(0), 
   fFirstPHOSParticle(-1), 
@@ -65,6 +67,28 @@ AliESD::~AliESD()
   fPmdTracks.Delete();
   fV0s.Delete();
   fCascades.Delete();
+  fKinks.Delete();
+  fV0MIs.Delete();
+}
+
+void AliESD::UpdateV0PIDs()
+{
+  //
+  //
+  //
+  Int_t nV0 = GetNumberOfV0MIs();
+  for (Int_t i=0;i<nV0;i++){
+    AliESDV0MI * v0 = GetV0MI(i);
+    AliESDtrack* tp = GetTrack(v0->fIndex[0]);
+    AliESDtrack* tm = GetTrack(v0->fIndex[1]);
+    if (!tm || !tp){
+      printf("BBBUUUUUUUGGGG\n");
+    }
+    Double_t pp[5],pm[5];
+    tp->GetESDpid(pp);
+    tm->GetESDpid(pm);
+    v0->UpdatePID(pp,pm);    
+  }
 }
 
 //______________________________________________________________________________
@@ -121,4 +145,6 @@ void AliESD::Print(Option_t *) const
   printf("                 pmd       %d\n", GetNumberOfPmdTracks());
   printf("                 v0        %d\n", GetNumberOfV0s());
   printf("                 cascades  %d\n)", GetNumberOfCascades());
+  printf("                 kinks     %d\n)", GetNumberOfKinks());
+  printf("                 V0MIs     %d\n)", GetNumberOfV0MIs());
 }
index d6352bab8579d4dc2b82cedb88fdbc745b70674a..7d59c331575b7c5aa59aded7f8d4540d7bcac8c9 100644 (file)
 #include "AliESDPmdTrack.h"
 #include "AliESDVertex.h"
 #include "AliESDcascade.h"
+#include "AliESDkink.h"
 #include "AliESDtrack.h"
 #include "AliESDHLTtrack.h"
 #include "AliESDv0.h"
+#include "AliESDV0MI.h"
 
 class AliESD : public TObject {
 public:
@@ -51,8 +53,11 @@ public:
     return (AliESDPmdTrack *)fPmdTracks.UncheckedAt(i);
   }
 
-  void AddTrack(const AliESDtrack *t) {
-    new(fTracks[fTracks.GetEntriesFast()]) AliESDtrack(*t);
+  Int_t  AddTrack(const AliESDtrack *t) {
+    AliESDtrack * track = new(fTracks[fTracks.GetEntriesFast()]) AliESDtrack(*t);
+    track->SetID(fTracks.GetEntriesFast()-1);
+    return  track->GetID();
+    
   }
   void AddHLTConfMapTrack(const AliESDHLTtrack *t) {
     new(fHLTConfMapTracks[fHLTConfMapTracks.GetEntriesFast()]) AliESDHLTtrack(*t);
@@ -73,6 +78,7 @@ public:
   void AddV0(const AliESDv0 *v) {
     new(fV0s[fV0s.GetEntriesFast()]) AliESDv0(*v);
   }
+  void UpdateV0PIDs();
 
   AliESDcascade *GetCascade(Int_t i) const {
     return (AliESDcascade *)fCascades.UncheckedAt(i);
@@ -81,6 +87,25 @@ public:
     new(fCascades[fCascades.GetEntriesFast()]) AliESDcascade(*c);
   }
 
+  AliESDkink *GetKink(Int_t i) const {
+    return (AliESDkink *)fKinks.UncheckedAt(i);
+  }
+  Int_t AddKink(const AliESDkink *c) {
+    AliESDkink * kink = new(fKinks[fKinks.GetEntriesFast()]) AliESDkink(*c);
+    kink->SetID(fKinks.GetEntriesFast());
+    return fKinks.GetEntriesFast()-1;
+  }
+
+  AliESDV0MI *GetV0MI(Int_t i) const {
+    return (AliESDV0MI *)fV0MIs.UncheckedAt(i);
+  }
+  Int_t AddV0MI(const AliESDV0MI *c) {
+    AliESDV0MI * v0 = new(fV0MIs[fV0MIs.GetEntriesFast()]) AliESDV0MI(*c);
+    v0->SetID(fV0MIs.GetEntriesFast()-1);
+    return fV0MIs.GetEntriesFast()-1;
+  }
+
+
   void SetVertex(AliESDVertex* vertex) {
     new(&fPrimaryVertex) AliESDVertex(*vertex);
   }
@@ -97,6 +122,8 @@ public:
   Int_t GetNumberOfPmdTracks() const {return fPmdTracks.GetEntriesFast();}
   Int_t GetNumberOfV0s()      const {return fV0s.GetEntriesFast();}
   Int_t GetNumberOfCascades() const {return fCascades.GetEntriesFast();}
+  Int_t GetNumberOfKinks() const {return fKinks.GetEntriesFast();}
+  Int_t GetNumberOfV0MIs() const {return fV0MIs.GetEntriesFast();}
   Int_t GetNumberOfPHOSParticles() const {return fPHOSParticles;}
   void  SetNumberOfPHOSParticles(Int_t part) { fPHOSParticles = part ; }
   void  SetFirstPHOSParticle(Int_t index) { fFirstPHOSParticle = index ; } 
@@ -145,6 +172,8 @@ protected:
   TClonesArray fPmdTracks;       // PMD ESD tracks
   TClonesArray fV0s;             // V0 vertices
   TClonesArray fCascades;        // Cascade vertices
+  TClonesArray fKinks;           // Kinks
+  TClonesArray fV0MIs;           // V0MI
   Int_t        fPHOSParticles;   // Number of PHOS particles (stored as fTracks)
   Int_t        fEMCALParticles;  // Number of EMCAL particles (stored as fTracks)
   Int_t        fFirstPHOSParticle; // First PHOS particle in the fTracks list 
index bc333e93fd065c679f18b1153f661ce321ba79bc..cda53c1364a1937d034139f1e783013c3b551e3c 100644 (file)
@@ -38,15 +38,16 @@ t2->Exec();
 
 //
 //some cuts definition
-TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.00005&&abs(MC.fVDist[2])<0.0005")
+TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01")
 //TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.5&&abs(MC.fVDist[2])<0.5")
-TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9");
+//TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9");
+TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<15");
 
 TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
 
 
 TCut crec("crec","fReconstructed==1");
-TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
+TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25");
 TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
 TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
 TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
@@ -54,20 +55,38 @@ TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fE
 
 AliESDComparisonDraw comp;  
 comp.SetIO(); 
+TFile f("genHits.root");
+TTree * treel = (TTree*)f.Get("HitLines");
+if (treel) comp->fTree->AddFriend(treel,"L");
 
 //
 //example
 comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
+comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
+comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
+comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
+comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
+
+comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
 
 TH1F his("his","his",100,0,20);
 TH1F hpools("hpools","hpools",100,-7,7);
 TH1F hfake("hfake","hfake",1000,0,150);
 TProfile profp0("profp0","profp0",20,-0.4,0.9)
 
-comp.DrawLogXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1&&abs(fPdg)==211"+cprim,"1",4,0.2,1.5,-0.06,0.06)
+comp.DrawLogXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06)
 comp.fRes->Draw();
 comp.fMean->Draw();  
-comp.Eff("fTPCinP0[3]","fRowsWithDigits>120&&abs(fPdg)==211"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
+
+comp.DrawLogXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06)
+comp.fRes->Draw();
+
+comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
+comp.fRes->Draw();
+
+comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio<0.1",10,0.2,1.5)
+comp.fRes->Draw();
+comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio>0.1",10,0.2,1.5)
 comp.fRes->Draw();
 
 comp.fTree->Draw("fESDTrack.fITSsignal/fESDTrack.fTPCsignal","fITSOn&&fTPCOn&&fESDTrack.fITSFakeRatio==0") 
@@ -75,15 +94,24 @@ comp.fTree->Draw("fESDTrack.fITSsignal/fESDTrack.fTPCsignal","fITSOn&&fTPCOn&&fE
 TH1F his("his","his",100,0,20);
 TH1F hpools("hpools","hpools",100,-7,7);
 
-TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,250); hdedx0->SetMarkerColor(2);
-TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,250); hdedx1->SetMarkerColor(3);
-TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,250); hdedx2->SetMarkerColor(4);
-TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,250); hdedx3->SetMarkerColor(5);
+TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(2);
+TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(3);
+TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(4);
+TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(6);
+
 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==211&&fITStrack.fN==6"+cprim) 
 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim) 
 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==321&&fITStrack.fN==6"+cprim) 
 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==11&&fITStrack.fN==6"+cprim) 
-hdedx0->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx3->Draw("same");
+
+
+comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&MC.fParticle.P()<2&&abs(fPdg)==211&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)") 
+comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&MC.fParticle.P()<2&&abs(fPdg)==2212&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)") 
+comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&MC.fParticle.P()<2&&abs(fPdg)==321&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)") 
+comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&MC.fParticle.P()<2&&abs(fPdg)==11&&fTRDtrack.fN>40&&abs(fESDTrack.fTRDLabel)==abs(fESDTrack.fTPCLabel)") 
+
+
+hdedx1->Draw(); hdedx0->Draw("same"); hdedx2->Draw("same"); hdedx3->Draw("same");
 
 comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
 
@@ -132,14 +160,84 @@ comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"
 #include "AliMagF.h"
 #include "AliESD.h"
 #include "AliESDtrack.h"
-#include "AliITStrackV2.h"
+#include "AliITStrackMI.h"
 #include "AliTRDtrack.h"
+#include "AliHelix.h"
+#include "AliESDVertex.h"
+#include "AliExternalTrackParam.h"
+#include "AliESDkink.h"
+#include "AliESDV0MI.h"
+
 #endif
 #include "AliGenInfo.h"
 #include "AliESDComparisonMI.h"
 
 
 
+
+void  AliESDRecInfo::UpdatePoints(AliESDtrack*track)
+{
+  //
+  //
+  Int_t iclusters[200];
+  Float_t density[160];
+  for (Int_t i=0;i<160;i++) density[i]=-1.;
+  fTPCPoints[0]= 160;
+  fTPCPoints[1] = -1;
+  //
+  if (fTPCPoints[0]<fTPCPoints[1]) return;
+  //  Int_t nclusters=track->GetTPCclusters(iclusters);
+
+  Int_t ngood=0;
+  Int_t undeff=0;
+  Int_t nall =0;
+  Int_t range=20;
+  for (Int_t i=0;i<160;i++){
+    Int_t last = i-range;
+    if (nall<range) nall++;
+    if (last>=0){
+      if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--;
+      if (iclusters[last]==-1) undeff--;
+    }
+    if (iclusters[i]>0&& (iclusters[i]&0x8000)==0)   ngood++;
+    if (iclusters[i]==-1) undeff++;
+    if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
+  }
+  Float_t maxdens=0;
+  Int_t indexmax =0;
+  for (Int_t i=0;i<160;i++){
+    if (density[i]<0) continue;
+    if (density[i]>maxdens){
+      maxdens=density[i];
+      indexmax=i;
+    }
+  }
+  //
+  //max dens point
+  fTPCPoints[3] = maxdens;
+  fTPCPoints[1] = indexmax;
+  //
+  // last point
+  for (Int_t i=indexmax;i<160;i++){
+    if (density[i]<0) continue;
+    if (density[i]<maxdens/2.) {
+      break;
+    }
+    fTPCPoints[2]=i;
+  }
+  //
+  // first point
+  for (Int_t i=indexmax;i>0;i--){
+    if (density[i]<0) continue;
+    if (density[i]<maxdens/2.) {
+      break;
+    }
+    fTPCPoints[0]=i;
+  }
+  //
+  if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1;
+}
+
 //
 //
 void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed)
@@ -149,6 +247,8 @@ void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstruct
   //calculates derived variables
   //  
   //
+  UpdatePoints(&fESDTrack);
+  fBestTOFmatch=1000;
   AliTrackReference * ref = &(info->fTrackRef);
   fTPCinR0[0] = info->fTrackRef.X();   
   fTPCinR0[1] = info->fTrackRef.Y();   
@@ -188,19 +288,74 @@ void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstruct
     fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
   }
   //
-  if (reconstructed==kFALSE){
-    fReconstructed = kFALSE;
-    fTPCOn = kFALSE;
-    fITSOn = kFALSE;
-    fTRDOn = kFALSE;
-    return;
-  }
+  for (Int_t i=0;i<4;i++) fStatus[i] =0;
+  fReconstructed = kFALSE;
+  fTPCOn = kFALSE;
+  fITSOn = kFALSE;
+  fTRDOn = kFALSE;  
+  if (reconstructed==kFALSE) return;
+
+  fLabels[0] = info->fLabel;
+  fLabels[1] = info->fPrimPart;
   fReconstructed = kTRUE;
   fTPCOn = ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
   fITSOn = ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
   fTRDOn = ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
-  
-  if (fTPCOn){
+  //
+  //  
+  if ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0){
+    fStatus[1] =3;
+  }
+  else{
+    if ((fESDTrack.GetStatus()&AliESDtrack::kTPCout)>0){
+      fStatus[1] =2;
+    }
+    else{
+      fStatus[1]=1;
+    }      
+  }
+  //
+  //
+  if ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0){
+    fStatus[2] =2;
+  }
+  else{
+    if ((fESDTrack.GetStatus()&AliESDtrack::kTRDout)>0){
+      fStatus[2] =1;
+    }
+  }
+  if ((fESDTrack.GetStatus()&AliESDtrack::kTRDStop)>0){
+    fStatus[2] =10;
+  }
+
+  //
+  //TOF 
+  // 
+  if (((fESDTrack.GetStatus()&AliESDtrack::kTOFout)>0)){
+    //
+    // best tof match
+    Double_t times[5];
+    fESDTrack.GetIntegratedTimes(times);    
+    for (Int_t i=0;i<5;i++){
+      if ( TMath::Abs(fESDTrack.GetTOFsignal()-times[i]) <TMath::Abs(fBestTOFmatch) ){
+       fBestTOFmatch = fESDTrack.GetTOFsignal()-times[i];
+      }
+    }
+    Int_t toflabel[3];
+    fESDTrack.GetTOFLabel(toflabel);
+    Bool_t toffake=kTRUE;
+    for (Int_t i=0;i<3;i++){
+      if (toflabel[i]<0) continue;
+      if (toflabel[i]== TMath::Abs(fESDTrack.GetLabel()))  toffake=kFALSE;     
+      fStatus[3]=1;
+    }
+    if (toffake) fStatus[3] =2;
+  }else{
+    fStatus[3]=0;
+  }
+
+
+  if (fStatus[1]>0 &&info->fNTPCRef>0){
     //TPC
     fESDTrack.GetInnerXYZ(fTPCinR1);
     fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
@@ -216,7 +371,8 @@ void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstruct
     }    
     Double_t cov[15], param[5],x;
     fESDTrack.GetInnerExternalCovariance(cov);
-    fESDTrack.GetInnerExternalParameters(x,param);    
+    fESDTrack.GetInnerExternalParameters(x,param);
+    //    if (x<50) return 
     //
     fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3];  //delta rfi
     fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
@@ -228,6 +384,7 @@ void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstruct
     fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
     fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
     Double_t sign = (param[4]>0)? 1.:-1; 
+    fSign =sign;
     fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(cov[14]);    
   }
   if (fITSOn){
@@ -264,6 +421,7 @@ void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstruct
     fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
     fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);    
     Double_t sign = (param[4]>0) ? 1:-1; 
+    fSign = sign;
     fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);    
   }
   
@@ -271,20 +429,261 @@ void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstruct
 
 
 void  AliESDRecV0Info::Update(Float_t vertex[3])
+{ 
+
+  if ( (fT1.fStatus[1]>0)&& (fT2.fStatus[1]>0)){
+    Float_t distance1,distance2;
+    Double_t xx[3],pp[3];
+    //
+    Double_t xd[3],pd[3],signd;
+    Double_t xm[3],pm[3],signm;
+    for (Int_t i=0;i<3;i++){
+      xd[i] = fT2.fTPCinR1[i];
+      pd[i] = fT2.fTPCinP1[i];
+      xm[i] = fT1.fTPCinR1[i];
+      pm[i] = fT1.fTPCinP1[i];
+    }
+    signd =  fT2.fSign<0 ? -1:1;
+    signm =  fT1.fSign<0 ? -1:1;
+
+    AliHelix dhelix1(xd,pd,signd);
+    dhelix1.GetMomentum(0,pp,0);
+    dhelix1.Evaluate(0,xx);      
+    // 
+    //  Double_t x2[3],p2[3];
+    //            
+    AliHelix mhelix(xm,pm,signm);    
+    //
+    //find intersection linear
+    //
+    Double_t phase[2][2],radius[2];
+    Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
+    Double_t delta1=10000,delta2=10000;  
+
+    if (points==1){
+      fMinR = TMath::Sqrt(radius[0]);
+    }
+    if (points==2){
+      fMinR =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
+    }
+    
+    if (points>0){
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    }
+    if (points==2){    
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    }
+    if (points==1){
+      fMinR = TMath::Sqrt(radius[0]);
+      fDistMinR = delta1;
+    }
+    if (points==2){
+      if (radius[0]<radius[1]){
+       fMinR = TMath::Sqrt(radius[0]);
+       fDistMinR = delta1;
+      }
+      else{
+       fMinR = TMath::Sqrt(radius[1]);
+       fDistMinR = delta2;
+      }
+    }
+    //
+    //
+    distance1 = TMath::Min(delta1,delta2);
+    //
+    //find intersection parabolic
+    //
+    points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+    delta1=10000,delta2=10000;  
+    
+    if (points>0){
+      dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    }
+    if (points==2){    
+      dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    }
+    
+    distance2 = TMath::Min(delta1,delta2);
+    if (delta1<delta2){
+      //get V0 info
+      dhelix1.Evaluate(phase[0][0],fXr);
+      dhelix1.GetMomentum(phase[0][0],fPdr);
+      mhelix.GetMomentum(phase[0][1],fPm);
+      dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
+      fRr = TMath::Sqrt(radius[0]);
+    }
+    else{
+      dhelix1.Evaluate(phase[1][0],fXr);
+      dhelix1.GetMomentum(phase[1][0], fPdr);
+      mhelix.GetMomentum(phase[1][1], fPm);
+      dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
+      fRr = TMath::Sqrt(radius[1]);
+    }
+    fDist1 = TMath::Sqrt(distance1);
+    fDist2 = TMath::Sqrt(distance2);      
+    
+    if (fDist2<10.5){
+      Double_t x,alpha,param[5],cov[15];
+      //
+      fT1.fESDTrack.GetInnerExternalParameters(x,param);
+      fT1.fESDTrack.GetInnerExternalCovariance(cov);
+      alpha = fT1.fESDTrack.GetInnerAlpha();
+      AliExternalTrackParam paramm(x,alpha,param,cov);
+      //
+      fT2.fESDTrack.GetInnerExternalParameters(x,param);
+      fT2.fESDTrack.GetInnerExternalCovariance(cov);
+      alpha = fT2.fESDTrack.GetInnerAlpha();
+      AliExternalTrackParam paramd(x,alpha,param,cov);
+    }    
+    //            
+    //   
+    
+    Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
+    Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
+    
+    Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
+    Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
+    vnorm2 = TMath::Sqrt(vnorm2);
+    Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
+    Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
+    pnorm2 = TMath::Sqrt(pnorm2);
+    
+    fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
+    fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
+    fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
+  }
+}
+
+////
+void  AliESDRecKinkInfo::Update()
 {
-  Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
-  Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
-
-  Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
-  Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
-  vnorm2 = TMath::Sqrt(vnorm2);
-  Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
-  Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
-  pnorm2 = TMath::Sqrt(pnorm2);
-  
-  fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
-  fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
-  fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
+
+  if ( (fT1.fTPCOn)&& (fT2.fTPCOn)){
+    //
+    // IF BOTH RECONSTRUCTED
+    Float_t distance1,distance2;
+    Double_t xx[3],pp[3];
+    //
+    Double_t xd[3],pd[3],signd;
+    Double_t xm[3],pm[3],signm;
+    for (Int_t i=0;i<3;i++){
+      xd[i] = fT2.fTPCinR1[i];
+      pd[i] = fT2.fTPCinP1[i];
+      xm[i] = fT1.fTPCinR1[i];
+      pm[i] = fT1.fTPCinP1[i];
+    }
+    signd =  fT2.fSign<0 ? -1:1;
+    signm =  fT1.fSign<0 ? -1:1;
+
+    AliHelix dhelix1(xd,pd,signd);
+    dhelix1.GetMomentum(0,pp,0);
+    dhelix1.Evaluate(0,xx);      
+    // 
+    //  Double_t x2[3],p2[3];
+    //            
+    AliHelix mhelix(xm,pm,signm);    
+    //
+    //find intersection linear
+    //
+    Double_t phase[2][2],radius[2];
+    Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
+    Double_t delta1=10000,delta2=10000;  
+
+    if (points==1){
+      fMinR = TMath::Sqrt(radius[0]);
+    }
+    if (points==2){
+      fMinR =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
+    }
+    
+    if (points>0){
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    }
+    if (points==2){    
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    }
+    if (points==1){
+      fMinR = TMath::Sqrt(radius[0]);
+      fDistMinR = delta1;
+    }
+    if (points==2){
+      if (radius[0]<radius[1]){
+       fMinR = TMath::Sqrt(radius[0]);
+       fDistMinR = delta1;
+      }
+      else{
+       fMinR = TMath::Sqrt(radius[1]);
+       fDistMinR = delta2;
+      }
+    }
+    //
+    //
+    distance1 = TMath::Min(delta1,delta2);
+    //
+    //find intersection parabolic
+    //
+    points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+    delta1=10000,delta2=10000;  
+    
+    if (points>0){
+      dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    }
+    if (points==2){    
+      dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    }
+    
+    distance2 = TMath::Min(delta1,delta2);
+    if (delta1<delta2){
+      //get V0 info
+      dhelix1.Evaluate(phase[0][0],fXr);
+      dhelix1.GetMomentum(phase[0][0],fPdr);
+      mhelix.GetMomentum(phase[0][1],fPm);
+      dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
+      fRr = TMath::Sqrt(radius[0]);
+    }
+    else{
+      dhelix1.Evaluate(phase[1][0],fXr);
+      dhelix1.GetMomentum(phase[1][0], fPdr);
+      mhelix.GetMomentum(phase[1][1], fPm);
+      dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
+      fRr = TMath::Sqrt(radius[1]);
+    }
+    fDist1 = TMath::Sqrt(distance1);
+    fDist2 = TMath::Sqrt(distance2);      
+    
+    if (fDist2<10.5){
+      Double_t x,alpha,param[5],cov[15];
+      //
+      fT1.fESDTrack.GetInnerExternalParameters(x,param);
+      fT1.fESDTrack.GetInnerExternalCovariance(cov);
+      alpha = fT1.fESDTrack.GetInnerAlpha();
+      AliExternalTrackParam paramm(x,alpha,param,cov);
+      //
+      fT2.fESDTrack.GetInnerExternalParameters(x,param);
+      fT2.fESDTrack.GetInnerExternalCovariance(cov);
+      alpha = fT2.fESDTrack.GetInnerAlpha();
+      AliExternalTrackParam paramd(x,alpha,param,cov);
+      /*
+      AliESDkink kink;
+      kink.Update(&paramm,&paramd);
+      //      kink.Dump();
+      Double_t diff  = kink.fRr-fRr;
+      Double_t diff2 = kink.fDist2-fDist2;
+      printf("Diff\t%f\t%f\n",diff,diff2);
+      */
+    }
+    
+    //            
+    //
+  }
 
 }
 
@@ -425,6 +824,8 @@ void ESDCmpTr::Reset()
   fEventNr = 0;
   fNEvents = 0;
   fTreeCmp = 0;
+  fTreeCmpKinks =0;
+  fTreeCmpV0 =0;
   //  fFnCmp = "cmpTracks.root";
   fFileGenTracks = 0;
   fDebug = 0;
@@ -451,6 +852,7 @@ Int_t ESDCmpTr::Exec()
     return 1;
    
   fNextTreeGenEntryToRead = 0;
+  fNextKinkToRead = 0;
   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
        eventNr++) {
@@ -458,17 +860,25 @@ Int_t ESDCmpTr::Exec()
     SetIO(fEventNr);
     fNParticles = gAlice->GetEvent(fEventNr);    
 
-    fIndexRecTracks = new Int_t[fNParticles*4];  //write at maximum 4 tracks corresponding to particle
+    fIndexRecTracks = new Int_t[fNParticles*20];  //write at maximum 4 tracks corresponding to particle
+    fIndexRecKinks  = new Int_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
+    fIndexRecV0  = new Int_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
+
     fFakeRecTracks = new Int_t[fNParticles];
     fMultiRecTracks = new Int_t[fNParticles];
-    for (Int_t i = 0; i<fNParticles; i++) {
-      fIndexRecTracks[4*i] = -1;
-      fIndexRecTracks[4*i+1] = -1;
-      fIndexRecTracks[4*i+2] = -1;
-      fIndexRecTracks[4*i+3] = -1;
+    fMultiRecKinks = new Int_t[fNParticles];
+    fMultiRecV0 = new Int_t[fNParticles];
 
+    for (Int_t i = 0; i<fNParticles; i++) {
+      for (Int_t j=0;j<20;j++){
+       fIndexRecTracks[20*i+j] = -1;
+       fIndexRecKinks[20*i+j]  = -1;
+       fIndexRecV0[20*i+j]  = -1;
+      }
       fFakeRecTracks[i] = 0;
       fMultiRecTracks[i] = 0;
+      fMultiRecKinks[i] = 0;
+      fMultiRecV0[i] = 0;      
     }
   
     cout<<"Start to process event "<<fEventNr<<endl;
@@ -478,11 +888,19 @@ Int_t ESDCmpTr::Exec()
 
     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
     if (TreeGenLoop(eventNr)>0) return 1;
+    BuildKinkInfo0(eventNr);
+    BuildV0Info(eventNr);
+    fRecArray->Delete();
+
     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
 
     delete [] fIndexRecTracks;
+    delete [] fIndexRecKinks;
+    delete [] fIndexRecV0;
     delete [] fFakeRecTracks;
     delete [] fMultiRecTracks;
+    delete [] fMultiRecKinks;
+    delete [] fMultiRecV0;
   }
 
   CloseOutputFile();
@@ -515,6 +933,29 @@ Bool_t ESDCmpTr::ConnectGenTree()
   fMCInfo = new AliMCInfo;
   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
   //
+  //
+  fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree");
+  if (!fTreeGenKinks) {
+    cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
+       <<fFnGenTracks<<endl;
+    //return kFALSE;
+  }
+  else{
+    fGenKinkInfo = new AliGenKinkInfo;
+    fTreeGenKinks->SetBranchAddress("MC",&fGenKinkInfo);
+  }
+
+  fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree");
+  if (!fTreeGenV0) {
+    cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
+       <<fFnGenTracks<<endl;
+    //return kFALSE;
+  }
+  else{
+    fGenV0Info = new AliGenV0Info;
+    fTreeGenV0->SetBranchAddress("MC",&fGenV0Info);
+  }
+  //
   if (fDebug > 1) {
     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
   }
@@ -530,24 +971,36 @@ void ESDCmpTr::CreateTreeCmp()
     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
     return;
   }
-
-
-  fTreeCmp    = new TTree("ESDcmpTracks","ESDcmpTracks");
   //
+  //
+  fTreeCmp    = new TTree("ESDcmpTracks","ESDcmpTracks");
   fMCInfo = new AliMCInfo;
   fRecInfo = new AliESDRecInfo;
-  //
   AliESDtrack * esdTrack = new AliESDtrack;
-  AliITStrackV2 * itsTrack = new AliITStrackV2;  
-  //
+  //  AliITStrackMI * itsTrack = new AliITStrackMI;  
   fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo);
   fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo);
   fTreeCmp->Branch("fESDTrack","AliESDtrack",&esdTrack);
-  //  fTreeCmp->Branch("ITS","AliITStrackV2",&itsTrack);
+  //  fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack);
   delete esdTrack;
   //
+  //
+  fTreeCmpKinks    = new TTree("ESDcmpKinks","ESDcmpKinks"); 
+  fGenKinkInfo     = new AliGenKinkInfo;
+  fRecKinkInfo     = new AliESDRecKinkInfo;
+  fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo);
+  fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo);
+  //
+  //
+  fTreeCmpV0       = new TTree("ESDcmpV0","ESDcmpV0"); 
+  fGenV0Info     = new AliGenV0Info;
+  fRecV0Info     = new AliESDRecV0Info;
+  fTreeCmpV0->Branch("MC.","AliGenV0Info",   &fGenV0Info);
+  fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info);
+  //
   fTreeCmp->AutoSave(); 
-
+  fTreeCmpKinks->AutoSave(); 
+  fTreeCmpV0->AutoSave(); 
 }
 ////////////////////////////////////////////////////////////////////////
 void ESDCmpTr::CloseOutputFile()  
@@ -582,23 +1035,35 @@ Int_t ESDCmpTr::TreeTLoop()
   //
   // loop over all ESD reconstructed tracks and store info in memory
   //
+  // + loop over all reconstructed kinks
   TStopwatch  timer;
   timer.Start();
   //  
   Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();  
+  Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks();
+  Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0MIs();
+  fSignedKinks = new Int_t[nKinks];
+  fSignedV0    = new Int_t[nV0MIs];
   //
-  fTracks      = new TObjArray(nEntries);
+  // load kinks to the memory
+  for (Int_t i=0; i<nKinks;i++){
+    AliESDkink * kink =fEvent->GetKink(i);
+    fSignedKinks[i]=0;
+    if (kink->fStatus<0) continue;
+  }
   //
-  //load tracks to the memory
-  for (Int_t i=0; i<nEntries;i++){
-    AliESDtrack * track =fEvent->GetTrack(i);
-    fTracks->AddAt(track,i);
+  for (Int_t i=0; i<nV0MIs;i++){
+    AliESDV0MI * v0MI =fEvent->GetV0MI(i);
+    fSignedV0[i]=0;
+    if (v0MI->fStatus<0) continue;
   }
+  
   //
   //
   AliESDtrack * track=0;
   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
-    track = (AliESDtrack*)fTracks->UncheckedAt(iEntry);
+    //track = (AliESDtrack*)fTracks->UncheckedAt(iEntry);
+    track = (AliESDtrack*)fEvent->GetTrack(iEntry);
     //
     Int_t label = track->GetLabel();
     Int_t absLabel = abs(label);
@@ -606,14 +1071,58 @@ Int_t ESDCmpTr::TreeTLoop()
       //      fIndexRecTracks[absLabel] =  iEntry;
       if (label < 0) fFakeRecTracks[absLabel]++;      
       if (fMultiRecTracks[absLabel]>0){
-       if (fMultiRecTracks[absLabel]<4)
-         fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] =  iEntry;      
+       if (fMultiRecTracks[absLabel]<20)
+         fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] =  iEntry;     
       }
       else      
-       fIndexRecTracks[absLabel*4] =  iEntry;
+       fIndexRecTracks[absLabel*20] =  iEntry;
       fMultiRecTracks[absLabel]++;
     }
+  }
+  // sort reconstructed kinks  
+  //
+  AliESDkink * kink=0;
+  for (Int_t iEntry=0; iEntry<nKinks;iEntry++){
+    kink = (AliESDkink*)fEvent->GetKink(iEntry);
+    if (!kink) continue;
+    //
+    Int_t label0 = TMath::Abs(kink->fLab[0]);
+    Int_t label1 = TMath::Abs(kink->fLab[1]);
+
+    Int_t absLabel = TMath::Min(label0,label1);
+    if (absLabel < fNParticles) {
+      if (fMultiRecKinks[absLabel]>0){
+       if (fMultiRecKinks[absLabel]<20)
+         fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] =  iEntry;       
+      }
+      else      
+       fIndexRecKinks[absLabel*20] =  iEntry;
+      fMultiRecKinks[absLabel]++;
+    }
   }  
+  // --sort reconstructed V0
+  //
+  AliESDV0MI * v0MI=0;
+  for (Int_t iEntry=0; iEntry<nV0MIs;iEntry++){
+    v0MI = fEvent->GetV0MI(iEntry);
+    if (!v0MI) continue;
+    //
+    Int_t label0 = TMath::Abs(v0MI->fLab[0]);
+    Int_t label1 = TMath::Abs(v0MI->fLab[1]);
+
+    Int_t absLabel = TMath::Min(label0,label1);
+    if (absLabel < fNParticles) {
+      if (fMultiRecV0[absLabel]>0){
+       if (fMultiRecV0[absLabel]<20)
+         fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] =  iEntry;     
+      }
+      else      
+       fIndexRecV0[absLabel*20] =  iEntry;
+      fMultiRecV0[absLabel]++;
+    }
+  }  
+
+
   printf("Time spended in TreeTLoop\n");
   timer.Print();
   
@@ -636,6 +1145,7 @@ Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
   TBranch * branch = fTreeCmp->GetBranch("RC");
   branch->SetAddress(&fRecInfo); // set all pointers
+  fRecArray = new TObjArray(fNParticles);
 
   while (entry < nParticlesTR) {
     fTreeGenTracks->GetEntry(entry);
@@ -647,6 +1157,7 @@ Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
     if (fDebug > 2 && fMCInfo->fLabel < 10) {
       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
     }
+    //    if (fMCInfo->fNTPCRef<1) continue; // not TPCref
     //
     fRecInfo->Reset();
     AliESDtrack * track=0;
@@ -654,29 +1165,60 @@ Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
     TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
     local.GetXYZ(fRecInfo->fTRLocalCoord);     
     //
-    if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
-      track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
+    if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) {
+      //track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
+      track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]);
       //
       //
       // find nearest track if multifound
-      if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
-       Double_t xyz[3];
-       track->GetInnerXYZ(xyz);
-       Float_t dz = TMath::Abs(local.Z()-xyz[2]);
-       for (Int_t i=1;i<4;i++){
-         if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
-           AliESDtrack * track2 = (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
-           track2->GetInnerXYZ(xyz);
-           Float_t dz2 = TMath::Abs(local.Z()-xyz[2]);
-           if  (TMath::Abs(dz2)<dz)
-             track = track2;              
+      //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge);
+      //
+      Int_t status = 0;
+      if  ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++;
+      if  ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++;
+      if  ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++;
+
+      //
+      if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){
+       //
+       Double_t p[3];
+       track->GetInnerPxPyPz(p);
+       Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
+       //
+       for (Int_t i=1;i<20;i++){
+         if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){
+           AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]);
+           if (!track2) continue;
+           //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; //           
+           //if (sign2<0) continue;
+           track2->GetInnerPxPyPz(p);
+           Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
+           /*
+           if (sign<0){
+             sign = sign2;
+             track = track2;
+             maxp = mom;
+             continue;
+           }
+           */
+           //
+           Int_t status2 = 0;
+           if  ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++;
+           if  ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++;
+           if  ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++;
+           if (status2<status) continue;
+           //
+           if (mom<maxp) continue;
+           maxp = mom;
+           track = track2;
+           //
          }
        }
       }        
       //
       fRecInfo->fESDTrack =*track; 
       if (track->GetITStrack())
-       fRecInfo->fITStrack = *((AliITStrackV2*)track->GetITStrack());
+       fRecInfo->fITStrack = *((AliITStrackMI*)track->GetITStrack());
       if (track->GetTRDtrack()){
        fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack());
       }
@@ -692,17 +1234,288 @@ Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
     else{
       fRecInfo->Update(fMCInfo,fParamTPC,kFALSE);
     }
-
+    fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel);
     fTreeCmp->Fill();
   }
   fTreeCmp->AutoSave();
-  fTracks->Delete();
+  //fTracks->Delete();
   printf("Time spended in TreeGenLoop\n");
   timer.Print();
   if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
 
   return 0;
 }
+
+
+
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+Int_t ESDCmpTr::BuildKinkInfo0(Int_t eventNr)
+{
+//
+// loop over all entries for a given event, find corresponding 
+// rec. track and store in the fTreeCmp
+//
+  TStopwatch timer;
+  timer.Start();
+  Int_t entry = fNextKinkToRead;
+  Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast();
+  cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" "
+      <<nParticlesTR<<" "<<fNextKinkToRead<<endl;
+  //
+  TBranch * branch = fTreeCmpKinks->GetBranch("RC.");
+  branch->SetAddress(&fRecKinkInfo); // set all pointers
+  
+  //
+  while (entry < nParticlesTR) {
+    fTreeGenKinks->GetEntry(entry);
+    entry++;
+    if (eventNr < fGenKinkInfo->fMCm.fEventNr) continue;
+    if (eventNr > fGenKinkInfo->fMCm.fEventNr) continue;;
+    //
+    fNextKinkToRead = entry-1;
+    //
+    //
+    AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCm.fLabel);
+    AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCd.fLabel);
+    fRecKinkInfo->fT1 = (*fRecInfo1);
+    fRecKinkInfo->fT2 = (*fRecInfo2);
+    fRecKinkInfo->fStatus =0;
+    if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1;
+    if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2;
+    if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1;
+    
+    if (fRecKinkInfo->fStatus==3){
+      fRecKinkInfo->Update();    
+    }
+    Int_t label =  TMath::Min(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
+    Int_t label2 = TMath::Max(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
+    
+    AliESDkink *kink=0;
+    fRecKinkInfo->fRecStatus   =0;
+    fRecKinkInfo->fMultiple    = fMultiRecKinks[label];
+    fRecKinkInfo->fKinkMultiple=0;
+    //
+    if (fMultiRecKinks[label]>0){
+
+      //      for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){
+      for (Int_t j=TMath::Min(fMultiRecKinks[label],20)-1;j>=0;j--){
+       Int_t index = fIndexRecKinks[label*20+j];
+       //AliESDkink *kink2  = (AliESDkink*)fKinks->At(index);
+       AliESDkink *kink2  = (AliESDkink*)fEvent->GetKink(index);
+       if (TMath::Abs(kink2->fLab[0])==label &&TMath::Abs(kink2->fLab[1])==label2) {
+         kink =kink2;
+         fRecKinkInfo->fKinkMultiple++;
+         fSignedKinks[index]=1;
+       }
+       if (TMath::Abs(kink2->fLab[1])==label &&TMath::Abs(kink2->fLab[0])==label2) {
+         kink =kink2;
+         fRecKinkInfo->fKinkMultiple++;
+         fSignedKinks[index]=1;
+       }
+      }
+    }
+    if (kink){
+      fRecKinkInfo->fKink = *kink;
+      fRecKinkInfo->fRecStatus=1;
+    }
+    fTreeCmpKinks->Fill();
+  }
+  //  Int_t nkinks = fKinks->GetEntriesFast();
+  Int_t nkinks = fEvent->GetNumberOfKinks();
+  for (Int_t i=0;i<nkinks;i++){
+    if (fSignedKinks[i]==0){
+      //      AliESDkink *kink  = (AliESDkink*)fKinks->At(i);
+      AliESDkink *kink  = (AliESDkink*)fEvent->GetKink(i);
+      if (!kink) continue;
+      //
+      fRecKinkInfo->fKink = *kink;
+      fRecKinkInfo->fRecStatus =-2;
+      //
+      AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[0]));
+      AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[1]));
+      if (fRecInfo1 && fRecInfo2){
+       fRecKinkInfo->fT1 = (*fRecInfo1);
+       fRecKinkInfo->fT2 = (*fRecInfo2);
+       fRecKinkInfo->fRecStatus =-1;
+      }
+      fTreeCmpKinks->Fill();
+    }
+  }
+
+
+  fTreeCmpKinks->AutoSave();
+  printf("Time spended in BuilKinkInfo Loop\n");
+  timer.Print();
+  if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl;
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+
+Int_t ESDCmpTr::BuildV0Info(Int_t eventNr)
+{
+//
+// loop over all entries for a given event, find corresponding 
+// rec. track and store in the fTreeCmp
+//
+  TStopwatch timer;
+  timer.Start();
+  Int_t entry = fNextV0ToRead;
+  Double_t nParticlesTR = fTreeGenV0->GetEntriesFast();
+  cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" "
+      <<nParticlesTR<<" "<<fNextV0ToRead<<endl;
+  //
+  TBranch * branch = fTreeCmpV0->GetBranch("RC.");
+  branch->SetAddress(&fRecV0Info); // set all pointers
+  const AliESDVertex * esdvertex = fEvent->GetVertex();
+  Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()};
+  
+  //
+  while (entry < nParticlesTR) {
+    fTreeGenV0->GetEntry(entry);
+    entry++;
+    if (eventNr < fGenV0Info->fMCm.fEventNr) continue;
+    if (eventNr > fGenV0Info->fMCm.fEventNr) continue;;
+    //
+    fNextV0ToRead = entry-1;
+    //
+    //
+    AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCm.fLabel);
+    AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCd.fLabel);
+    if (fGenV0Info->fMCm.fCharge*fGenV0Info->fMCd.fCharge>0) continue;  // interactions
+    if (!fRecInfo1 || !fRecInfo2) continue;
+    fRecV0Info->fT1 = (*fRecInfo1);
+    fRecV0Info->fT2 = (*fRecInfo2);
+    fRecV0Info->fV0Status =0;
+    if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1;
+    if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2;
+
+    if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1;
+
+
+    if (abs(fRecV0Info->fV0Status)==3){
+      fRecV0Info->Update(vertex);
+      {
+       //
+       // TPC V0 Info
+       Double_t x,alpha, param[5],cov[15];
+       fRecV0Info->fT1.fESDTrack.GetInnerExternalParameters(x,param);
+       fRecV0Info->fT1.fESDTrack.GetInnerExternalCovariance(cov);
+       alpha = fRecV0Info->fT1.fESDTrack.GetInnerAlpha();
+       AliExternalTrackParam paramP(x,alpha,param,cov);
+       //
+       fRecV0Info->fT2.fESDTrack.GetInnerExternalParameters(x,param);
+       fRecV0Info->fT2.fESDTrack.GetInnerExternalCovariance(cov);
+       alpha = fRecV0Info->fT2.fESDTrack.GetInnerAlpha();
+       AliExternalTrackParam paramM(x,alpha,param,cov);
+       //
+       fRecV0Info->fV0tpc.SetM(paramM);
+       fRecV0Info->fV0tpc.SetP(paramP);
+       Double_t pid1[5],pid2[5];
+       fRecV0Info->fT1.fESDTrack.GetESDpid(pid1);
+       fRecV0Info->fT1.fESDTrack.GetESDpid(pid2);
+       //
+       fRecV0Info->fV0tpc.UpdatePID(pid1,pid2);
+       fRecV0Info->fV0tpc.Update(vertex);
+       //
+       //
+       fRecV0Info->fT1.fESDTrack.GetExternalParameters(x,param);
+       fRecV0Info->fT1.fESDTrack.GetExternalCovariance(cov);
+       alpha = fRecV0Info->fT1.fESDTrack.GetAlpha();
+       new (&paramP) AliExternalTrackParam(x,alpha,param,cov);
+       //
+       fRecV0Info->fT2.fESDTrack.GetExternalParameters(x,param);
+       fRecV0Info->fT2.fESDTrack.GetExternalCovariance(cov);
+       alpha = fRecV0Info->fT2.fESDTrack.GetAlpha();
+       new (&paramM) AliExternalTrackParam(x,alpha,param,cov);
+       //
+       fRecV0Info->fV0its.SetM(paramM);
+       fRecV0Info->fV0its.SetP(paramP);
+       fRecV0Info->fV0its.UpdatePID(pid1,pid2);
+       fRecV0Info->fV0its.Update(vertex);
+
+      }
+      if (TMath::Abs(fGenV0Info->fMCm.fPdg)==11 &&TMath::Abs(fGenV0Info->fMCd.fPdg)==11){
+       if (fRecV0Info->fDist2>10){
+         fRecV0Info->Update(vertex);
+       }
+       if (fRecV0Info->fDist2>10){
+         fRecV0Info->Update(vertex);
+       }
+      }
+    }   
+    //
+    // take the V0 from reconstruction
+    Int_t label =  TMath::Min(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);
+    Int_t label2 = TMath::Max(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);    
+    AliESDV0MI *v0MI=0;
+    fRecV0Info->fRecStatus   =0;
+    fRecV0Info->fMultiple    = fMultiRecV0[label];
+    fRecV0Info->fV0Multiple=0;
+    //
+    if (fMultiRecV0[label]>0){
+
+      //      for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){
+      for (Int_t j=TMath::Min(fMultiRecV0[label],20)-1;j>=0;j--){
+       Int_t index = fIndexRecV0[label*20+j];
+       AliESDV0MI *v0MI2  = fEvent->GetV0MI(index);
+       if (TMath::Abs(v0MI2->fLab[0])==label &&TMath::Abs(v0MI2->fLab[1])==label2) {
+         v0MI =v0MI2;
+         fRecV0Info->fV0Multiple++;
+         fSignedV0[index]=1;
+       }
+       if (TMath::Abs(v0MI2->fLab[1])==label &&TMath::Abs(v0MI2->fLab[0])==label2) {
+         v0MI =v0MI2;
+         fRecV0Info->fV0Multiple++;
+         fSignedV0[index]=1;
+       }
+      }
+    }
+    if (v0MI){
+      fRecV0Info->fV0rec = *v0MI;
+      fRecV0Info->fRecStatus=1;
+    }
+
+    fTreeCmpV0->Fill();
+  }
+  //
+  // write fake v0s
+
+  Int_t nV0MIs = fEvent->GetNumberOfV0MIs();
+  for (Int_t i=0;i<nV0MIs;i++){
+    if (fSignedV0[i]==0){
+      AliESDV0MI *v0MI  = fEvent->GetV0MI(i);
+      if (!v0MI) continue;
+      //
+      fRecV0Info->fV0rec = *v0MI;
+      fRecV0Info->fV0Status  =-10;
+      fRecV0Info->fRecStatus =-2;
+      //
+      AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[0]));
+      AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[1]));
+      if (fRecInfo1 && fRecInfo2){
+       fRecV0Info->fT1 = (*fRecInfo1);
+       fRecV0Info->fT2 = (*fRecInfo2);
+       fRecV0Info->fRecStatus =-1;
+      }
+      fTreeCmpV0->Fill();
+    }
+  }
+
+
+
+  fTreeCmpV0->AutoSave();
+  fRecArray->Delete();
+  printf("Time spended in BuilV0Info Loop\n");
+  timer.Print();
+  if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<endl;
+  return 0;
+}
 ////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////
 
index 10899a28d920518e28da8055ce3e91a123ff4079..9f8ff5a297b5f3f1409fe99c897fbc19eeef5da0 100644 (file)
@@ -14,8 +14,9 @@ class AliESDRecInfo: public TObject {
 public:
   AliESDRecInfo(){}
   ~AliESDRecInfo(){}
-  //
+  void UpdatePoints(AliESDtrack* track);
   void Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed);
+  Float_t    fTPCPoints[4]; //start , biggest end points + max density
   Double_t fTPCinR0[5];   //generated position of the track at inner tpc - radius [3] and fi [4]
   Double_t fTPCinR1[5];   //reconstructed postion of the track           - radius [3] and fi [
   Double_t fTPCinP0[5];   //generated position of the track at inner tpc
@@ -33,16 +34,20 @@ public:
   Double_t fITSDelta[5];  // deltas
   Double_t fITSPools[5];  // pools
   AliESDtrack fESDTrack;          // tpc track
-  AliITStrackV2 fITStrack;        //its track
+  AliITStrackMI fITStrack;        //its track
   AliTRDtrack fTRDtrack;        //its track
+  Float_t fBestTOFmatch;        //best matching between times
   Float_t fTRLocalCoord[3];       //local coordinates of the track ref.
   Int_t   fReconstructed;         //flag if track was reconstructed
   Int_t fFake;             // fake track
   Int_t fMultiple;         // number of reconstructions
   Bool_t fTPCOn;           // TPC refitted inward
+  Int_t  fStatus[4];        // status -0 not found - 1 -only in - 2 -in-out -3 -in -out-refit
   Bool_t fITSOn;           // ITS refitted inward
   Bool_t fTRDOn;           // ITS refitted inward
   Float_t fDeltaP;          //delta of momenta
+  Double_t fSign;
+  Int_t fLabels[2];
   void Reset();
   ClassDef(AliESDRecInfo,1)  // container for 
 };
@@ -73,6 +78,8 @@ public:
   Double_t       fPdr[3];    //momentum at vertex daughter  - according approx at DCA
   Double_t       fXr[3];     //rec. position according helix
   //
+  Double_t       fMinR;     // minimum radius in rphi intersection
+  Double_t       fDistMinR; // distance at minimal radius
   Double_t       fPm[3];    //momentum at the vertex mother
   Double_t       fAngle[3]; //three angles
   Double_t       fRr;       // rec position of the vertex 
@@ -80,12 +87,51 @@ public:
   Float_t        fPointAngleFi; //point angle fi
   Float_t        fPointAngleTh; //point angle theta
   Float_t        fPointAngle;   //point angle full
+  Int_t          fV0Status;       // status of the kink
+  AliESDV0MI     fV0tpc;           // Vo information from reconsturction according TPC
+  AliESDV0MI     fV0its;           // Vo information from reconsturction according ITS
+  AliESDV0MI     fV0rec;           // V0 information form the reconstruction
+  Int_t          fMultiple;
+  Int_t          fV0Multiple;
+  Int_t          fRecStatus;    // status form the reconstuction
   ClassDef(AliESDRecV0Info,1)   // container for  
 };
 
 ClassImp(AliESDRecV0Info)
 
 
+class AliESDRecKinkInfo: public TObject {
+public:
+  void Update();
+  AliESDRecInfo  fT1;      //track1
+  AliESDRecInfo  fT2;      //track2  
+  AliESDkink     fKink;    //kink
+  Double_t       fDist1;    //info about closest distance according closest MC - linear DCA
+  Double_t       fDist2;    //info about closest distance parabolic DCA
+  Double_t       fInvMass;  //reconstructed invariant mass -
+  //
+  Double_t       fPdr[3];    //momentum at vertex daughter  - according approx at DCA
+  Double_t       fXr[3];     //rec. position according helix
+  //
+  Double_t       fPm[3];    //momentum at the vertex mother
+  Double_t       fAngle[3]; //three angles
+  Double_t       fRr;       // rec position of the vertex 
+  Double_t       fMinR;     // minimum radius in rphi intersection
+  Double_t       fDistMinR; // distance at minimal radius
+  Int_t          fLab[2];   //MC label of the partecle
+  Float_t        fPointAngleFi; //point angle fi
+  Float_t        fPointAngleTh; //point angle theta
+  Float_t        fPointAngle;   //point angle full
+  Int_t          fStatus;       //status -tracks 
+  Int_t          fRecStatus;    //kink -status- 0 - not found  1-good -  fake
+  Int_t          fMultiple;
+  Int_t          fKinkMultiple;
+  ClassDef(AliESDRecKinkInfo,1)   // container for  
+};
+
+ClassImp(AliESDRecKinkInfo)
+
+
 
 ////////////////////////////////////////////////////////////////////////
 // 
@@ -112,6 +158,8 @@ public:
   Bool_t ConnectGenTree();
   Int_t TreeGenLoop(Int_t eventNr);
   Int_t TreeTLoop();
+  Int_t BuildKinkInfo0(Int_t eventNr); // build kink info 0
+  Int_t BuildV0Info(Int_t eventNr); // build kink info 0
   void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
   void SetNEvents(Int_t i) {fNEvents = i;}
   void SetDebug(Int_t level) {fDebug = level;}
@@ -129,10 +177,14 @@ private:
   char  fFnCmp[1000];                   //! output file name with cmp tracks
   TFile *fFileCmp;                //! output file with cmp tracks
   TTree *fTreeCmp;                //! output tree with cmp tracks
+  TTree *fTreeCmpKinks;                //! output tree with cmp Kinks
+  TTree *fTreeCmpV0;                //! output tree with cmp V0
   //
   char  fFnGenTracks[1000];             //! input file name with gen tracks
   TFile *fFileGenTracks;
   TTree *fTreeGenTracks;
+  TTree *fTreeGenKinks;            // tree with gen kinks
+  TTree *fTreeGenV0;            // tree with gen V0
   //
   //
   Int_t  fDirection;
@@ -144,17 +196,30 @@ private:
   Int_t *fFakeRecTracks;          //! number of fake tracks
   Int_t *fMultiRecTracks;         //! number of multiple reconstructions
   //
-  TObjArray * fTracks;            //!container with tracks 
+  Int_t *fIndexRecKinks;         //! index of particle label in treeesd
+  Int_t *fMultiRecKinks;         //! number of multiple reconstructions
+  Int_t *fSignedKinks;                //! indicator that kink was not fake
+  //
+  Int_t *fIndexRecV0;         //! index of particle label in treeesd
+  Int_t *fMultiRecV0;         //! number of multiple reconstructions
+  Int_t *fSignedV0;                //! indicator that kink was not fake
+  //
+  TObjArray *fRecArray;           // container with rec infos
   AliESD *fEvent;                 //!event
-
   //
   AliTPCParam* fParamTPC;         //! AliTPCParam
   Int_t fNParticles;              //! number of particles in the input tree genTracks
   Int_t fDebug;                   //! debug flag  
   Int_t fNextTreeGenEntryToRead;    //! last entry already read from genTracks tree
+  Int_t fNextKinkToRead;            //! last entry already read from genKinks tree
+  Int_t fNextV0ToRead;            //! last entry already read from genV0 tree
   //
   AliMCInfo*  fMCInfo;           //! MC information writen per particle
+  AliGenKinkInfo* fGenKinkInfo;      //! MC information writen per Kink
+  AliGenV0Info* fGenV0Info;      //! MC information writen per Kink
   AliESDRecInfo*  fRecInfo;          //! Rec. information writen per particle
+  AliESDRecKinkInfo* fRecKinkInfo;    //! reconstructed kink info
+  AliESDRecV0Info* fRecV0Info;    //! reconstructed kink info
   //
 
   ClassDef(ESDCmpTr,1)    // class which creates and fills tree with ESDGenTrack objects
diff --git a/STEER/AliESDV0MI.cxx b/STEER/AliESDV0MI.cxx
new file mode 100644 (file)
index 0000000..3e9326c
--- /dev/null
@@ -0,0 +1,216 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+//    Origin: Marian Ivanov marian.ivanov@cern.ch
+//-------------------------------------------------------------------------
+
+#include <Riostream.h>
+#include <TMath.h>
+#include <TPDGCode.h>
+#include "AliESDV0MI.h"
+#include "AliHelix.h"
+
+
+ClassImp(AliESDV0MI)
+
+AliESDV0MI::AliESDV0MI(){
+  //
+  //Dafault constructor
+  //
+  fID =0;
+  fDist1  =-1;
+  fDist2  =-1;
+  fRr     =-1;
+  fStatus = 0;
+  fRow0   =-1;  
+}
+
+void AliESDV0MI::SetP(const AliExternalTrackParam & paramp)  {
+  //
+  // set mother
+  //
+  fParamP   = paramp;
+}
+
+void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
+  //
+  //set daughter
+  //
+  fParamM = paramm;
+
+}
+  
+void  AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5])
+{
+  //
+  // set PID hypothesy
+  //
+  // norm PID to 1
+  Float_t sump =0;
+  Float_t summ =0;
+  for (Int_t i=0;i<5;i++){
+    fRP[i]=pidp[i];
+    sump+=fRP[i];
+    fRM[i]=pidm[i];
+    summ+=fRM[i];
+  }
+  for (Int_t i=0;i<5;i++){
+    fRP[i]/=sump;
+    fRM[i]/=summ;
+  }
+}
+
+Float_t AliESDV0MI::GetProb(UInt_t p1, UInt_t p2){
+  //
+  //
+  //
+  //
+  return TMath::Max(fRP[p1]+fRM[p2], fRP[p2]+fRM[p1]);
+}
+
+Float_t AliESDV0MI::GetEffMass(UInt_t p1, UInt_t p2){
+  //
+  // calculate effective mass
+  //
+  const Float_t pmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
+                     4.93599999999999983e-01, 9.38270000000000048e-01};
+  if (p1>4) return -1;
+  if (p2>4) return -1;
+  Float_t mass1 = pmass[p1]; 
+  Float_t mass2 = pmass[p2];   
+  Double_t *m1 = fPP;
+  Double_t *m2 = fPM;
+  //
+  if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
+    m1 = fPM;
+    m2 = fPP;
+  }
+  //
+  Float_t e1    = TMath::Sqrt(mass1*mass1+
+                              m1[0]*m1[0]+
+                              m1[1]*m1[1]+
+                              m1[2]*m1[2]);
+  Float_t e2    = TMath::Sqrt(mass2*mass2+
+                              m2[0]*m2[0]+
+                              m2[1]*m2[1]+
+                              m2[2]*m2[2]);  
+  Float_t mass =  
+    (m2[0]+m1[0])*(m2[0]+m1[0])+
+    (m2[1]+m1[1])*(m2[1]+m1[1])+
+    (m2[2]+m1[2])*(m2[2]+m1[2]);
+  
+  mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
+  return mass;
+}
+
+void  AliESDV0MI::Update(Float_t vertex[3])
+{
+  //
+  // updates Kink Info
+  //
+  Float_t distance1,distance2;
+  //
+  AliHelix phelix(fParamP);
+  AliHelix mhelix(fParamM);    
+  //
+  //find intersection linear
+  //
+  Double_t phase[2][2],radius[2];
+  Int_t  points = phelix.GetRPHIintersections(mhelix, phase, radius,200);
+  Double_t delta1=10000,delta2=10000;  
+  
+  if (points>0){
+    phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+  }
+  if (points==2){    
+    phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+  }
+  distance1 = TMath::Min(delta1,delta2);
+  //
+  //find intersection parabolic
+  //
+  points = phelix.GetRPHIintersections(mhelix, phase, radius);
+  delta1=10000,delta2=10000;  
+  Double_t d1=1000.,d2=10000.;
+  if (points>0){
+    phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    phelix.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    Double_t xd[3],xm[3];
+    phelix.Evaluate(phase[0][0],xd);
+    mhelix.Evaluate(phase[0][1],xm);
+    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]);
+  }
+  if (points==2){    
+    phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    phelix.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    Double_t xd[3],xm[3];
+    phelix.Evaluate(phase[1][0],xd);
+    mhelix.Evaluate(phase[1][1],xm);
+    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]);
+  }
+  //
+  distance2 = TMath::Min(delta1,delta2);
+  if (delta1<delta2){
+    //get V0 info
+    Double_t xd[3],xm[3];
+    phelix.Evaluate(phase[0][0],xd);
+    mhelix.Evaluate(phase[0][1], xm);
+    fXr[0] = 0.5*(xd[0]+xm[0]);
+    fXr[1] = 0.5*(xd[1]+xm[1]);
+    fXr[2] = 0.5*(xd[2]+xm[2]);
+    //
+    phelix.GetMomentum(phase[0][0],fPP);
+    mhelix.GetMomentum(phase[0][1],fPM);
+    phelix.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
+    fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
+  }
+  else{
+    Double_t xd[3],xm[3];
+    phelix.Evaluate(phase[1][0],xd);
+    mhelix.Evaluate(phase[1][1], xm);
+    fXr[0] = 0.5*(xd[0]+xm[0]);
+    fXr[1] = 0.5*(xd[1]+xm[1]);
+    fXr[2] = 0.5*(xd[2]+xm[2]);
+    //
+    phelix.GetMomentum(phase[1][0], fPP);
+    mhelix.GetMomentum(phase[1][1], fPM);
+    phelix.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
+    fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
+  }
+  fDist1 = TMath::Sqrt(TMath::Min(d1,d2));
+  fDist2 = TMath::Sqrt(distance2);      
+  //            
+  //
+  Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
+  Float_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
+  Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
+  Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
+  vnorm2 = TMath::Sqrt(vnorm2);
+  Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
+  Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
+  pnorm2 = TMath::Sqrt(pnorm2);  
+  fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
+  fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
+  fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
+  //
+}
+
diff --git a/STEER/AliESDV0MI.h b/STEER/AliESDV0MI.h
new file mode 100644 (file)
index 0000000..c7d731a
--- /dev/null
@@ -0,0 +1,111 @@
+#ifndef ALIESDV0MI_H
+#define ALIESDV0MI_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+//                          ESD V0 Vertex Class
+//          This class is part of the Event Summary Data set of classes
+//    Origin: Marian Ivanov marian.ivanov@cern.ch
+//-------------------------------------------------------------------------
+
+#include <AliESDv0.h>
+#include "AliExternalTrackParam.h"
+#include <TPDGCode.h>
+
+class AliESDtrack;
+
+class AliESDV0MI :  public AliESDv0 {
+public:
+  AliESDV0MI();             //constructor
+  //
+  void SetP(const AliExternalTrackParam & paramp); 
+  void SetM(const AliExternalTrackParam & paramd);
+  void UpdatePID(Double_t pidp[5], Double_t pidm[5]);
+  Float_t GetEffMass(UInt_t p1, UInt_t p2);
+  Float_t GetProb(UInt_t p1, UInt_t p2);
+  void Update(Float_t vertex[3]);            //update
+  void SetID(Int_t id){fID =id;}
+  Int_t GetID(){ return fID;}
+  AliExternalTrackParam fParamP;
+  AliExternalTrackParam fParamM;
+  Float_t        fRP[5];         // combined pid positive
+  Float_t        fRM[5];         // combined pid positive
+  Int_t          fID;
+  Int_t          fLab[2];     //MC label of the partecle
+  Int_t          fIndex[2];   //reconstructed labels of the tracks
+  //
+  //  
+  Double_t       fDist1;    //info about closest distance according closest MC - linear DCA
+  Double_t       fDist2;    //info about closest distance parabolic DCA
+  //
+  Double_t       fPP[3];    //momentum  positive   - according approx at DCA
+  Double_t       fPM[3];    //momentum negative
+  //
+  Double_t       fXr[3];      //rec. position according helix
+  Double_t       fAngle[3];   //three angles
+  Double_t       fRr;         //rec position of the vertex 
+  Int_t          fStatus;       //status 
+  Int_t          fRow0;         // critical layer
+  Int_t          fOrder[3]; //order of the vertex 
+  //  quality information
+  Double_t       fDistNorm; //normalized  DCA
+  Double_t       fDistSigma; //sigma of distance
+  Float_t        fChi2Before;   //chi2 of the tracks before V0
+  Float_t        fNBefore;      // number of possible points before V0
+  Float_t        fChi2After;   // chi2 of the tracks after V0
+  Float_t        fNAfter;      // number of possible points after V0
+  Float_t        fPointAngleFi; //point angle fi
+  Float_t        fPointAngleTh; //point angle theta
+  Float_t        fPointAngle;   //point angle full
+
+
+
+  ClassDef(AliESDV0MI,1)      // ESD V0 vertex
+};
+
+
+
+
+#endif
+
+/*
+
+class AliITSRecV0Info: public AliV0vertex {
+  friend class AliITStrackerMI;
+protected:
+  AliITSRecV0Info();
+  AliITSRecV0Info(const AliITStrackV2 &neg, const AliITStrackV2 &pos);
+  void Update(Float_t vertex[3], Float_t mass1, Float_t mass2);
+  Double_t       fDist1;    //info about closest distance according closest MC - linear DCA
+  Double_t       fDist2;    //info about closest distance parabolic DCA
+  Double_t       fDistNorm; //normalized  DCA
+  Double_t       fDistSigma; //sigma of distance
+  Double_t       fInvMass;  //reconstructed invariant mass -
+  //
+  Double_t       fPdr[3];    //momentum at vertex daughter  - according approx at DCA
+  Double_t       fXr[3];     //rec. position according helix
+  //
+  Double_t       fPm[3];    //momentum at the vertex mother
+  Double_t       fAngle[3]; //three angles
+  Double_t       fRr;       // rec position of the vertex 
+  Int_t          fLab[2];   //MC label of the particle
+  Int_t          fIndex[2]; //indexes of the tracks
+  Int_t          fOrder[3]; //order of the vertex 
+  Float_t        fChi2Before;   //chi2 of the tracks before V0
+  Float_t        fNBefore;      // number of possible points before V0
+  //
+  Float_t        fChi2After;   // chi2 of the tracks after V0
+  Float_t        fNAfter;      // number of possible points after V0
+  //
+  Float_t        fPointAngleFi; //point angle fi
+  Float_t        fPointAngleTh; //point angle theta
+  Float_t        fPointAngle;   //point angle full
+  ClassDef(AliITSRecV0Info,1)  // container for  
+};
+
+*/
diff --git a/STEER/AliESDkink.cxx b/STEER/AliESDkink.cxx
new file mode 100644 (file)
index 0000000..e478873
--- /dev/null
@@ -0,0 +1,161 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+//    Origin: Marian Ivanov marian.ivanov@cern.ch
+//-------------------------------------------------------------------------
+
+#include <Riostream.h>
+#include <TMath.h>
+#include <TPDGCode.h>
+#include "AliESDkink.h"
+#include "AliHelix.h"
+
+
+ClassImp(AliESDkink)
+
+AliESDkink::AliESDkink(){
+  //
+  //Dafault constructor
+  //
+  fDist1  =-1;
+  fDist2  =-1;
+  fRr     =-1;
+  fStatus = 0;
+  fRow0   =-1;
+  fTPCdensity2[0][0]=-1;
+  fTPCdensity2[0][1]=-1;
+
+}
+
+void AliESDkink::SetMother(const AliExternalTrackParam & pmother)  {
+  //
+  // set mother
+  //
+  fParamMother   = pmother;
+}
+
+void AliESDkink::SetDaughter(const AliExternalTrackParam & pdaughter){
+  //
+  //set daughter
+  //
+  fParamDaughter = pdaughter;
+
+}
+  
+void  AliESDkink::Update()
+{
+  //
+  // updates Kink Info
+  //
+  Float_t distance1,distance2;
+  //
+  AliHelix dhelix1(fParamDaughter);
+  AliHelix mhelix(fParamMother);    
+  //
+  //find intersection linear
+  //
+  Double_t phase[2][2],radius[2];
+  Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
+  Double_t delta1=10000,delta2=10000;  
+  
+  if (points>0){
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+  }
+  if (points==2){    
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+  }
+  distance1 = TMath::Min(delta1,delta2);
+  //
+  //find intersection parabolic
+  //
+  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+  delta1=10000,delta2=10000;  
+  Double_t d1=1000.,d2=10000.;
+  if (points>0){
+    dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    Double_t xd[3],xm[3];
+    dhelix1.Evaluate(phase[0][0],xd);
+    mhelix.Evaluate(phase[0][1],xm);
+    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]);
+  }
+  if (points==2){    
+    dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    Double_t xd[3],xm[3];
+    dhelix1.Evaluate(phase[1][0],xd);
+    mhelix.Evaluate(phase[1][1],xm);
+    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]);
+  }
+  //
+  distance2 = TMath::Min(delta1,delta2);
+  if (delta1<delta2){
+    //get V0 info
+    //    dhelix1.Evaluate(phase[0][0],fXr);
+    Double_t xd[3],xm[3];
+    dhelix1.Evaluate(phase[0][0],xd);
+    mhelix.Evaluate(phase[0][1], xm);
+    fXr[0] = 0.5*(xd[0]+xm[0]);
+    fXr[1] = 0.5*(xd[1]+xm[1]);
+    fXr[2] = 0.5*(xd[2]+xm[2]);
+    //
+    dhelix1.GetMomentum(phase[0][0],fPdr);
+    mhelix.GetMomentum(phase[0][1],fPm);
+    dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
+    //fRr = TMath::Sqrt(radius[0]);
+    fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
+  }
+  else{
+    //dhelix1.Evaluate(phase[1][0],fXr);
+    Double_t xd[3],xm[3];
+    dhelix1.Evaluate(phase[1][0],xd);
+    mhelix.Evaluate(phase[1][1], xm);
+    fXr[0] = 0.5*(xd[0]+xm[0]);
+    fXr[1] = 0.5*(xd[1]+xm[1]);
+    fXr[2] = 0.5*(xd[2]+xm[2]);
+    //
+    dhelix1.GetMomentum(phase[1][0], fPdr);
+    mhelix.GetMomentum(phase[1][1], fPm);
+    dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
+    //    fRr = TMath::Sqrt(radius[1]); 
+    fRr = TMath::Sqrt(fXr[0]*fXr[0]+fXr[1]*fXr[1]);
+  }
+  fDist1 = TMath::Sqrt(TMath::Min(d1,d2));
+  fDist2 = TMath::Sqrt(distance2);      
+  //            
+  //
+
+}
+
+Float_t AliESDkink::GetTPCDensityFactor() const
+{
+  //
+  //
+  return fTPCdensity[0][0]+fTPCdensity[1][1]-TMath::Max(fTPCdensity[0][1],Float_t(0.0))-TMath::Max(fTPCdensity[1][0],Float_t(0.0)); 
+}
+
+Float_t AliESDkink::GetQt() const
+{
+  Float_t dmomentum = TMath::Sqrt(fPdr[0]*fPdr[0]+fPdr[1]*fPdr[1]+fPdr[2]*fPdr[2]);
+  return TMath::Sin(fAngle[2])*dmomentum;
+}
diff --git a/STEER/AliESDkink.h b/STEER/AliESDkink.h
new file mode 100644 (file)
index 0000000..99cf28a
--- /dev/null
@@ -0,0 +1,60 @@
+#ifndef ALIESDKINK_H
+#define ALIESDKINK_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+//                          ESD V0 Vertex Class
+//          This class is part of the Event Summary Data set of classes
+//    Origin: Marian Ivanov marian.ivanov@cern.ch
+//-------------------------------------------------------------------------
+
+#include <TObject.h>
+#include "AliExternalTrackParam.h"
+#include <TPDGCode.h>
+
+class AliESDtrack;
+
+class AliESDkink : public TObject {
+public:
+  AliESDkink();             //constructor
+  //
+  void SetID(Int_t id){fID=id;}
+  Int_t GetID(){return fID;}
+  void SetMother(const AliExternalTrackParam & pmother); 
+  void SetDaughter(const AliExternalTrackParam & pdaughter);
+  void Update();            //update
+  Float_t GetTPCDensityFactor() const;
+  Float_t GetQt() const;    
+  //  
+  Int_t          fID;       // kink ID
+  AliExternalTrackParam fParamDaughter;
+  AliExternalTrackParam fParamMother;
+  Double_t       fDist1;    //info about closest distance according closest MC - linear DCA
+  Double_t       fDist2;    //info about closest distance parabolic DCA
+  //
+  Double_t       fPdr[3];    //momentum at vertex daughter  - according approx at DCA
+  Double_t       fXr[3];     //rec. position according helix
+  //
+  Double_t       fPm[3];    //momentum at the vertex mother
+  Double_t       fAngle[3]; //three angles
+  Double_t       fRr;       // rec position of the vertex 
+  Int_t          fLab[2];   //MC label of the partecle
+  Int_t          fIndex[2]; //reconstructed labels of the tracks
+  Int_t          fStatus;       //status 
+  Float_t        fTPCdensity[2][2];  //tpc cluster density before and after kink
+  Float_t        fTPCdensity2[2][2];  //tpc cluster density before and after kink - after second iteration
+  Float_t        fShapeFactor;       // tpc clusters shape factor
+  Int_t          fRow0;              // critical pad row number
+  Int_t          fMultiple[2];       //how many times the track's were used
+  Float_t        fZm[2];                // z at the middle of chamber
+  Float_t        fFi[2];
+  ClassDef(AliESDkink,1)      // ESD V0 vertex
+};
+
+#endif
+
+
index 9510f11024ee0208ac5a22098e4c119cf44a6259..3b7d4108c5cf3af5a9fbc99bb5733db824b55b0f 100644 (file)
@@ -66,6 +66,7 @@ fRICHsignal(-1)
   //
   // The default ESD constructor 
   //
+  fID =0;
   for (Int_t i=0; i<kSPECIES; i++) {
     fTrackTime[i]=0.;
     fR[i]=1.;
@@ -93,7 +94,12 @@ fRICHsignal(-1)
   }
   for (i=0; i<6; i++)  { fITSindex[i]=0; }
   for (i=0; i<180; i++){ fTPCindex[i]=0; }
+  for (i=0; i<3;i++)   { fKinkIndexes[i]=0;}
+  for (i=0; i<3;i++)   { fV0Indexes[i]=-1;}
   for (i=0; i<130; i++) { fTRDindex[i]=0; }
+  for (Int_t i=0;i<4;i++) {fTPCPoints[i]=-1;}
+  for (Int_t i=0;i<3;i++) {fTOFLabel[i]=-1;}
+  for (Int_t i=0;i<10;i++) {fTOFInfo[i]=-1;}
   fTPCLabel = 0;
   fTRDLabel = 0;
   fITSLabel = 0;
@@ -107,6 +113,7 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):TObject(track){
   //
   //copy constructor
   //
+  fID = track.fID;
   fFlags = track.fFlags;
   fLabel =track.fLabel;
   fTrackLength =track.fTrackLength;
@@ -162,6 +169,9 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):TObject(track){
   fTPCsignal=track.fTPCsignal;      
   for (Int_t i=0;i<kSPECIES;i++) fTPCr[i]=track.fTPCr[i]; 
   fTPCLabel=track.fTPCLabel;       
+  for (Int_t i=0;i<4;i++) {fTPCPoints[i]=track.fTPCPoints[i];}
+  for (Int_t i=0; i<3;i++)   { fKinkIndexes[i]=track.fKinkIndexes[i];}
+  for (Int_t i=0; i<3;i++)   { fV0Indexes[i]=track.fV0Indexes[i];}
   //
   fTRDchi2=track.fTRDchi2;        
   fTRDncls=track.fTRDncls;       
@@ -176,6 +186,8 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):TObject(track){
   fTOFindex=track.fTOFindex;       
   fTOFsignal=track.fTOFsignal;      
   for (Int_t i=0;i<kSPECIES;i++) fTOFr[i]=track.fTOFr[i];
+  for (Int_t i=0;i<3;i++) fTOFLabel[i]=track.fTOFLabel[i];
+  for (Int_t i=0;i<10;i++) fTOFInfo[i]=track.fTOFInfo[i];
   //
   for (Int_t i=0;i<3;i++) fPHOSpos[i]=track.fPHOSpos[i]; 
   fPHOSsignal=track.fPHOSsignal; 
@@ -258,7 +270,7 @@ Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags) {
   case kTPCin: case kTPCrefit:
     fTPCLabel = t->GetLabel();
     fIalpha=fRalpha;
-    fIx=fRx;
+    fIx=fRx;    
     {
       Int_t i;
       for (i=0; i<5; i++) fIp[i]=fRp[i];
@@ -339,8 +351,7 @@ Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags) {
        t->GetExternalCovariance(fXc); //can be done better
     }
   case kTRDin: case kTRDrefit:
-    fTRDLabel = t->GetLabel();
+    fTRDLabel = t->GetLabel(); 
     fTRDncls=t->GetNumberOfClusters();
     fTRDchi2=t->GetChi2();
     for (Int_t i=0;i<fTRDncls;i++) fTRDindex[i]=t->GetClusterIndex(i);
@@ -767,12 +778,36 @@ void AliESDtrack::SetTOFpid(const Double_t *p) {
   SetStatus(AliESDtrack::kTOFpid);
 }
 
+//_______________________________________________________________________
+void AliESDtrack::SetTOFLabel(const Int_t *p) {  
+  // Sets  (in TOF)
+  for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
+}
+
 //_______________________________________________________________________
 void AliESDtrack::GetTOFpid(Double_t *p) const {
   // Gets probabilities of each particle type (in TOF)
   for (Int_t i=0; i<kSPECIES; i++) p[i]=fTOFr[i];
 }
 
+//_______________________________________________________________________
+void AliESDtrack::GetTOFLabel(Int_t *p) const {
+  // Gets (in TOF)
+  for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
+}
+
+//_______________________________________________________________________
+void AliESDtrack::GetTOFInfo(Float_t *info) const {
+  // Gets (in TOF)
+  for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
+}
+
+//_______________________________________________________________________
+void AliESDtrack::SetTOFInfo(Float_t*info) {
+  // Gets (in TOF)
+  for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
+}
+
 
 
 //_______________________________________________________________________
index a7ff877e14dcc3a80570d6af966e1ce6cb0cd755..9f235a22b625f21e9964f8201a9d28d61d1ef6d3 100644 (file)
@@ -20,7 +20,9 @@ class AliESDtrack : public TObject {
 public:
   AliESDtrack();
   AliESDtrack(const AliESDtrack& track);
-  virtual ~AliESDtrack();  
+  virtual ~AliESDtrack();
+  void SetID(Int_t id) { fID =id;}
+  Int_t GetID(){ return fID;}
   void SetStatus(ULong_t flags) {fFlags|=flags;}
   void ResetStatus(ULong_t flags) {fFlags&=~flags;}
   Bool_t UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags);
@@ -79,10 +81,15 @@ public:
 
   void SetTPCpid(const Double_t *p);
   void GetTPCpid(Double_t *p) const;
+  void SetTPCPoints(Float_t points[4]){for (Int_t i=0;i<4;i++) fTPCPoints[i]=points[i];}
+  void SetKinkIndexes(Int_t points[3]) {for (Int_t i=0;i<3;i++) fKinkIndexes[i] = points[i];}
+  void SetV0Indexes(Int_t points[3]) {for (Int_t i=0;i<3;i++) fV0Indexes[i] = points[i];}
   Float_t GetTPCsignal() const {return fTPCsignal;}
   Float_t GetTPCchi2() const {return fTPCchi2;}
   Int_t GetTPCclusters(Int_t *idx) const;
   Int_t GetTPCLabel() const {return fTPCLabel;}
+  Int_t GetKinkIndex(Int_t i) const { return fKinkIndexes[i];}
+  Int_t GetV0Index(Int_t i) const { return fV0Indexes[i];}
   const TBits& GetTPCClusterMap() const {return fTPCClusterMap;}
   
   void SetTRDpid(const Double_t *p);
@@ -91,6 +98,7 @@ public:
   Float_t GetTRDsignal() const {return fTRDsignal;}
   Float_t GetTRDchi2() const {return fTRDchi2;}
   Int_t GetTRDclusters(UInt_t *idx) const;
+  Int_t GetTRDncls() const {return fTRDncls;}
   void    SetTRDpid(Int_t iSpecies, Float_t p);
   Float_t GetTRDpid(Int_t iSpecies) const;
   Int_t GetTRDLabel() const {return fTRDLabel;}
@@ -101,7 +109,11 @@ public:
   Float_t GetTOFsignal() const {return fTOFsignal;}
   Float_t GetTOFchi2() const {return fTOFchi2;}
   void    SetTOFpid(const Double_t *p);
+  void    SetTOFLabel(const Int_t *p);
   void    GetTOFpid(Double_t *p) const;
+  void    GetTOFLabel(Int_t *p) const;
+  void    GetTOFInfo(Float_t *info) const;
+  void    SetTOFInfo(Float_t *info);
   UInt_t  GetTOFcluster() const {return fTOFindex;}
   void  SetTOFcluster(UInt_t index) {fTOFindex=index;}
   
@@ -159,7 +171,7 @@ public:
 protected:
   ULong_t   fFlags;        // Reconstruction status flags 
   Int_t     fLabel;        // Track label
-
+  Int_t     fID;           // Unique ID of the track
   Float_t   fTrackLength;         // Track length
   Float_t   fTrackTime[kSPECIES]; // TOFs estimated by the tracking
   Float_t   fR[kSPECIES];         // combined "detector response probability"
@@ -216,11 +228,14 @@ protected:
   // TPC related track information
   Float_t fTPCchi2;        // chi2 in the TPC
   Int_t   fTPCncls;        // number of clusters assigned in the TPC
-  UInt_t  fTPCindex[180];  //! indices of the assigned TPC clusters
+  Int_t  fTPCindex[180];  //! indices of the assigned TPC clusters
   TBits   fTPCClusterMap;  // Map of clusters, one bit per padrow; 1 if has a cluster on given padrow
   Float_t fTPCsignal;      // detector's PID signal
   Float_t fTPCr[kSPECIES]; // "detector response probabilities" (for the PID)
   Int_t   fTPCLabel;       // label according TPC
+  Float_t fTPCPoints[4];   // TPC points -first, max. dens, last and max density
+  Int_t   fKinkIndexes[3]; // array of indexes of posible kink candidates 
+  Int_t   fV0Indexes[3]; // array of indexes of posible kink candidates 
   // TRD related track information
   Float_t fTRDchi2;        // chi2 in the TRD
   Int_t   fTRDncls;        // number of clusters assigned in the TRD
@@ -235,9 +250,10 @@ protected:
   UInt_t  fTOFindex;       // index of the assigned TOF cluster
   Float_t fTOFsignal;      // detector's PID signal
   Float_t fTOFr[kSPECIES]; // "detector response probabilities" (for the PID)
-
+  Int_t   fTOFLabel[3];       // TOF label 
+  Float_t fTOFInfo[10];       //! TOF informations
   // PHOS related track information 
-  Float_t fPHOSpos[3]; //position localised by PHOS in global coordinate system
+  Float_t fPHOSpos[3]; // position localised by PHOS in global coordinate system
   Float_t fPHOSsignal; // energy measured by PHOS
   Float_t fPHOSr[kSPECIESN]; // PID information from PHOS
 
@@ -250,7 +266,7 @@ protected:
   Float_t fRICHsignal;     // detector's PID signal (beta for RICH)
   Float_t fRICHr[kSPECIES];// "detector response probabilities" (for the PID)
        
-  ClassDef(AliESDtrack,8)  //ESDtrack 
+  ClassDef(AliESDtrack,9)  //ESDtrack 
 };
 
 #endif 
diff --git a/STEER/AliExternalTrackParam.cxx b/STEER/AliExternalTrackParam.cxx
new file mode 100644 (file)
index 0000000..3c7c73f
--- /dev/null
@@ -0,0 +1,262 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+// track parameters in "external" format                                     //
+//                                                                           //
+// The track parameters are:                                                //
+// - local y coordinate                                                      //
+// - local z coordinate                                                      //
+// - sin of azimutal angle                                                   //
+// - tan of dip angle                                                        //
+// - charge/pt                                                               //
+// The parametrisation is given at the local x coordinate fX and the         //
+// azimuthal angle fAlpha.                                                   //
+//                                                                           //
+// The external parametrisation can be used to exchange track parameters     //
+// between different detectors.                                              //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+
+#include "AliExternalTrackParam.h"
+#include "AliKalmanTrack.h"
+
+ClassImp(AliExternalTrackParam)
+
+
+//_____________________________________________________________________________
+AliExternalTrackParam::AliExternalTrackParam()
+{
+// default constructor
+
+  fX = fAlpha = 0.;
+  for (Int_t i = 0; i < 5; i++) fParam[i] = 0;
+  for (Int_t i = 0; i < 15; i++) fCovar[i] = 0;
+}
+
+//_____________________________________________________________________________
+AliExternalTrackParam::AliExternalTrackParam(Double_t x, Double_t alpha, 
+                                            const Double_t param[5], 
+                                            const Double_t covar[15])
+{
+// create external track parameters from given arguments
+
+  fX = x;
+  fAlpha = alpha;
+  for (Int_t i = 0; i < 5; i++) fParam[i] = param[i];
+  for (Int_t i = 0; i < 15; i++) fCovar[i] = covar[i];
+}
+
+AliExternalTrackParam::AliExternalTrackParam(const AliKalmanTrack& track)
+{
+  //
+  //
+  track.GetExternalParameters(fX,fParam);
+  track.GetExternalCovariance(fCovar);
+  fAlpha = track.GetAlpha();
+}
+
+
+//_____________________________________________________________________________
+const Double_t* AliExternalTrackParam::GetParameter() const
+{
+// get a pointer to the array of track parameters
+
+  return fParam;
+}
+
+//_____________________________________________________________________________
+const Double_t* AliExternalTrackParam::GetCovariance() const
+{
+// get a pointer to the array of the track parameter covariance matrix
+
+  return fCovar;
+}
+
+//_____________________________________________________________________________
+AliExternalTrackParam* AliExternalTrackParam::CreateExternalParam() const
+{
+// copy this instance
+
+  return new AliExternalTrackParam(fX, fAlpha, fParam, fCovar);
+}
+
+//_____________________________________________________________________________
+void AliExternalTrackParam::ResetCovariance(Double_t factor,
+                                           Bool_t clearOffDiagonal)
+{
+// reset the covariance matrix ("forget" track history)
+
+  Int_t k = 0;
+  for (Int_t i = 0; i < 5; i++) {
+    for (Int_t j = 0; j < i; j++) {  // off diagonal elements
+      if (clearOffDiagonal) {
+       fCovar[k++] = 0;
+      } else {
+       fCovar[k++] *= factor;
+      }
+    }
+    fCovar[k++] *= factor;     // diagonal elements
+  }
+}
+
+
+//_____________________________________________________________________________
+Bool_t AliExternalTrackParam::PropagateTo(Double_t /*x*/, Double_t* /*length*/)
+{
+// Propagate the track parameters to the given x coordinate assuming vacuum.
+// If length is not NULL, the change of track length is added to it.
+//
+// NOT IMPLEMENTED for this class
+
+  return kFALSE;
+}
+
+//_____________________________________________________________________________
+Bool_t AliExternalTrackParam::RotateTo(Double_t /*alpha*/)
+{
+// Rotate the reference axis for the parametrisation to the given angle.
+//
+// NOT IMPLEMENTED for this class
+
+  return kFALSE;
+}
+
+//_____________________________________________________________________________
+Bool_t AliExternalTrackParam::CorrectForMaterial(Double_t /*dAngle*/, 
+                                                Double_t /*dPrel*/)
+{
+// Take into account material effects assuming:
+//   dAngle2: mean multiple scattering angle in rad squared
+//   dPrel  : mean relative momentum gain (if > 0) or loss (if < 0)
+//
+// NOT IMPLEMENTED for this class
+
+  return kFALSE;
+}
+
+//_____________________________________________________________________________
+Bool_t AliExternalTrackParam::GetProlongationAt(Double_t /*x*/, 
+                                               Double_t& /*y*/, 
+                                               Double_t& /*z*/) const
+{
+// Get the local y and z coordinates at the given x value
+//
+// NOT IMPLEMENTED for this class
+
+  return kFALSE;
+}
+
+//_____________________________________________________________________________
+Double_t AliExternalTrackParam::GetXAtVertex(Double_t /*x*/, 
+                                            Double_t /*y*/) const
+{
+// Get the x coordinate at the given vertex (x,y)
+//
+// NOT IMPLEMENTED for this class
+
+  return 0;
+}
+
+
+//_____________________________________________________________________________
+Double_t AliExternalTrackParam::GetPredictedChi2(const AliCluster* /*cluster*/)
+{
+// calculate the chi2 contribution of the given cluster
+//
+// NOT IMPLEMENTED for this class
+
+  return -1;
+}
+
+//_____________________________________________________________________________
+Bool_t AliExternalTrackParam::Update(const AliCluster* /*cluster*/)
+{
+// update the track parameters using the position and error 
+// of the given cluster
+//
+// NOT IMPLEMENTED for this class
+
+  return kFALSE;
+}
+
+
+//_____________________________________________________________________________
+Double_t AliExternalTrackParam::SigmaPhi() const
+{
+// get the error of the azimuthal angle
+
+  return TMath::Sqrt(TMath::Abs(fCovar[5] / (1. - fParam[2]*fParam[2])));
+}
+
+//_____________________________________________________________________________
+Double_t AliExternalTrackParam::SigmaTheta() const
+{
+// get the error of the polar angle
+
+  return TMath::Sqrt(TMath::Abs(fCovar[9])) / (1. + fParam[3]*fParam[3]);
+}
+
+//_____________________________________________________________________________
+Double_t AliExternalTrackParam::SigmaPt() const
+{
+// get the error of the transversal component of the momentum
+
+  return TMath::Sqrt(fCovar[14]) / TMath::Abs(fParam[4]);
+}
+
+//_____________________________________________________________________________
+TVector3 AliExternalTrackParam::Momentum() const
+{
+// get the momentum vector
+
+  Double_t phi = TMath::ASin(fParam[2]) + fAlpha;
+  Double_t pt = 1. / TMath::Abs(fParam[4]);
+  return TVector3(pt * TMath::Cos(phi), 
+                 pt * TMath::Sin(phi), 
+                 pt * fParam[3]);
+}
+
+//_____________________________________________________________________________
+TVector3 AliExternalTrackParam::Position() const
+{
+// get the current spatial position in global coordinates
+
+  return TVector3(fX * TMath::Cos(fAlpha) - fParam[0] * TMath::Sin(fAlpha),
+                 fX * TMath::Sin(fAlpha) + fParam[0] * TMath::Cos(fAlpha),
+                 fParam[1]);
+}
+
+
+//_____________________________________________________________________________
+void AliExternalTrackParam::Print(Option_t* /*option*/) const
+{
+// print the parameters and the covariance matrix
+
+  printf("AliExternalTrackParam: x = %-12g  alpha = %-12g\n", fX, fAlpha);
+  printf("  parameters: %12g %12g %12g %12g %12g\n",
+        fParam[0], fParam[1], fParam[2], fParam[3], fParam[4]);
+  printf("  covariance: %12g\n", fCovar[0]);
+  printf("              %12g %12g\n", fCovar[1], fCovar[2]);
+  printf("              %12g %12g %12g\n", fCovar[3], fCovar[4], fCovar[5]);
+  printf("              %12g %12g %12g %12g\n", 
+        fCovar[6], fCovar[7], fCovar[8], fCovar[9]);
+  printf("              %12g %12g %12g %12g %12g\n", 
+        fCovar[10], fCovar[11], fCovar[12], fCovar[13], fCovar[14]);
+}
diff --git a/STEER/AliExternalTrackParam.h b/STEER/AliExternalTrackParam.h
new file mode 100644 (file)
index 0000000..1120571
--- /dev/null
@@ -0,0 +1,56 @@
+#ifndef ALIEXTERNALTRACKPARAM_H
+#define ALIEXTERNALTRACKPARAM_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+
+#include "AliTrackParam.h"
+class AliKalmanTrack;
+
+class AliExternalTrackParam: public AliTrackParam {
+ public:
+  AliExternalTrackParam();
+  AliExternalTrackParam(Double_t x, Double_t alpha, 
+                       const Double_t param[5], const Double_t covar[15]);
+  AliExternalTrackParam(const AliKalmanTrack& track);
+
+  virtual const Double_t* GetParameter() const;
+  virtual const Double_t* GetCovariance() const;
+  virtual Double_t     X() const {return fX;};
+  virtual Double_t     Alpha() const {return fAlpha;};
+  virtual AliExternalTrackParam* CreateExternalParam() const;
+  virtual void         ResetCovariance(Double_t factor = 10.,
+                                      Bool_t clearOffDiagonal = kTRUE);
+  virtual Double_t     Y() const {return fParam[0];};
+  virtual Double_t     Z() const {return fParam[1];};
+
+  virtual Bool_t       PropagateTo(Double_t x, Double_t* length = NULL);
+  virtual Bool_t       RotateTo(Double_t alpha);
+  virtual Bool_t       CorrectForMaterial(Double_t dAngle2, Double_t dPrel);
+  virtual Bool_t       GetProlongationAt(Double_t x, Double_t& y, 
+                                        Double_t& z) const;
+  virtual Double_t     GetXAtVertex(Double_t x = 0, Double_t y = 0) const;
+
+  virtual Double_t     GetPredictedChi2(const AliCluster* cluster);
+  virtual Bool_t       Update(const AliCluster* cluster);
+
+  virtual Double_t     SigmaPhi() const;
+  virtual Double_t     SigmaTheta() const;
+  virtual Double_t     SigmaPt() const;
+  virtual TVector3     Momentum() const;
+  virtual TVector3     Position() const;
+
+  virtual void         Print(Option_t* option = "") const;
+
+ private:
+  Double_t             fX;          // x coordinate for the parametrisation
+  Double_t             fAlpha;      // azimuthal angle for the parametrisation
+  Double_t             fParam[5];   // track parameter (y, z, sin(azimuthal angel), tan(dip angle), 1/pt)
+  Double_t             fCovar[15];  // track parameter covariance
+
+  ClassDef(AliExternalTrackParam, 1)
+};
+
+#endif
index ab21613989fd265dd3c9a3311c6fb4c4415020a6..a26684b22c88520107c50274daa1a09b5a6dbc03 100644 (file)
@@ -49,6 +49,9 @@ t->Exec();
 #include "TCanvas.h"
 #include "TPad.h"
 #include "TF1.h"
+#include "TView.h"
+#include "TGeometry.h"
+#include "TPolyLine3D.h"
 
 //ALIROOT includes
 #include "AliRun.h"
@@ -62,10 +65,13 @@ t->Exec();
 #include "AliTPCParamSR.h"
 #include "AliTracker.h"
 #include "AliMagF.h"
+#include "AliHelix.h"
+#include "AliPoints.h"
+
 #endif
 #include "AliGenInfo.h" 
 //
-//
+// 
 
 AliTPCParam * GetTPCParam(){
   AliTPCParamSR * par = new AliTPCParamSR;
@@ -98,6 +104,92 @@ Float_t TPCBetheBloch(Float_t bg)
   return ((Float_t)((kp2-aa-bb)*kp1/aa));
 }
 
+AliPointsMI::AliPointsMI(){
+  fN=0;
+  fX=0;
+  fY=0;
+  fZ=0;  
+  fCapacity = 0;
+  fLabel0=0;
+  fLabel1=0;
+}
+
+
+AliPointsMI::AliPointsMI(Int_t n, Float_t *x,Float_t *y, Float_t *z){
+  //
+  //
+  fN=n;
+  fCapacity = 1000+2*n;
+  fX= new Float_t[n];
+  fY= new Float_t[n];
+  fZ= new Float_t[n];
+  memcpy(fX,x,n*sizeof(Float_t));
+  memcpy(fY,y,n*sizeof(Float_t));
+  memcpy(fZ,z,n*sizeof(Float_t));
+  fLabel0=0;
+  fLabel1=0;
+}
+
+void AliPointsMI::Reset()
+{
+  fN=0;
+}
+
+void AliPointsMI::Reset(AliDetector * det, Int_t particle){
+  //
+  // write points from detector points
+  //
+  Reset();
+  if (!det) return;
+  TObjArray *array = det->Points();
+  if (!array) return;
+  for (Int_t i=0;i<array->GetEntriesFast();i++){
+    AliPoints * points = (AliPoints*) array->At(i);
+    if (!points) continue;
+    if (points->GetIndex()!=particle) continue;
+    Int_t npoints = points->GetN();
+    if (npoints<2) continue;
+    Int_t delta = npoints/100;
+    if (delta<1) delta=1;
+    if (delta>10) delta=10;
+    Int_t mypoints = npoints/delta;
+    //
+    fN = mypoints;
+    if (fN>fCapacity){
+      fCapacity = 1000+2*fN;
+      delete []fX;
+      delete []fY;
+      delete []fZ;
+      fX = new Float_t[fCapacity];
+      fY = new Float_t[fCapacity];
+      fZ = new Float_t[fCapacity]; 
+    }   
+    Float_t *xyz = points->GetP();
+    for (Int_t ipoint=0;ipoint<mypoints;ipoint++){
+      Int_t index = 3*ipoint*delta;
+      fX[ipoint]=0;
+      fY[ipoint]=0;
+      fZ[ipoint]=0;
+      if (index+2<npoints*3){
+       fX[ipoint] = xyz[index];
+       fY[ipoint] = xyz[index+1];
+       fZ[ipoint] = xyz[index+2];
+      }
+    }
+  }
+  fLabel0 = particle;
+}
+
+
+AliPointsMI::~AliPointsMI(){
+  fN=0;
+  fCapacity =0;
+  delete[] fX;
+  delete[]fY;
+  delete []fZ;
+  fX=0;fY=0;fZ=0;  
+}
+
 
 
 
@@ -108,6 +200,7 @@ AliMCInfo::AliMCInfo()
   fTPCReferences  = new TClonesArray("AliTrackReference",10);
   fITSReferences  = new TClonesArray("AliTrackReference",10);
   fTRDReferences  = new TClonesArray("AliTrackReference",10);
+  fTOFReferences  = new TClonesArray("AliTrackReference",10);
   fTRdecay.SetTrack(-1);
 }
 
@@ -122,6 +215,10 @@ AliMCInfo::~AliMCInfo()
   if (fTRDReferences){    
     delete fTRDReferences;  
   }
+  if (fTOFReferences){    
+    delete fTOFReferences;  
+  }
+
 }
 
 
@@ -139,7 +236,9 @@ void AliMCInfo::Update()
   //Float_t rlast=0;
   fNTPCRef = fTPCReferences->GetEntriesFast();
   fNITSRef = fITSReferences->GetEntriesFast();
-
+  fNTRDRef = fTRDReferences->GetEntriesFast();
+  fNTOFRef = fTOFReferences->GetEntriesFast();
+  
   for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
     AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
     //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
@@ -199,25 +298,294 @@ void AliMCInfo::Update()
 }
 
 /////////////////////////////////////////////////////////////////////////////////
+/*
 void AliGenV0Info::Update()
 {
   fMCPd[0] = fMCd.fParticle.Px();
   fMCPd[1] = fMCd.fParticle.Py();
   fMCPd[2] = fMCd.fParticle.Pz();
   fMCPd[3] = fMCd.fParticle.P();
+  //
   fMCX[0]  = fMCd.fParticle.Vx();
   fMCX[1]  = fMCd.fParticle.Vy();
   fMCX[2]  = fMCd.fParticle.Vz();
   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
+  //
   fPdg[0]    = fMCd.fParticle.GetPdgCode();
   fPdg[1]    = fMCm.fParticle.GetPdgCode();
   //
   fLab[0]    = fMCd.fParticle.GetUniqueID();
   fLab[1]    = fMCm.fParticle.GetUniqueID();
+  //  
+}
+*/
 
+void AliGenV0Info::Update(Float_t vertex[3])
+{
+  fMCPd[0] = fMCd.fParticle.Px();
+  fMCPd[1] = fMCd.fParticle.Py();
+  fMCPd[2] = fMCd.fParticle.Pz();
+  fMCPd[3] = fMCd.fParticle.P();
+  //
+  fMCX[0]  = fMCd.fParticle.Vx();
+  fMCX[1]  = fMCd.fParticle.Vy();
+  fMCX[2]  = fMCd.fParticle.Vz();
+  fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
+  //
+  fPdg[0]    = fMCd.fParticle.GetPdgCode();
+  fPdg[1]    = fMCm.fParticle.GetPdgCode();
+  //
+  fLab[0]    = fMCd.fParticle.GetUniqueID();
+  fLab[1]    = fMCm.fParticle.GetUniqueID();
+  //
+  //
+  //
+  Double_t x1[3],p1[3];
+  TParticle & pdaughter = fMCd.fParticle;
+  x1[0] = pdaughter.Vx();      
+  x1[1] = pdaughter.Vy();
+  x1[2] = pdaughter.Vz();
+  p1[0] = pdaughter.Px();
+  p1[1] = pdaughter.Py();
+  p1[2] = pdaughter.Pz();
+  Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
+  AliHelix dhelix1(x1,p1,sign);
+  //
+  //
+  Double_t x2[3],p2[3];            
+  //
+  TParticle & pmother = fMCm.fParticle;
+  x2[0] = pmother.Vx();      
+  x2[1] = pmother.Vy();      
+  x2[2] = pmother.Vz();      
+  p2[0] = pmother.Px();
+  p2[1] = pmother.Py();
+  p2[2] = pmother.Pz();
+  //
+  //
+  sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
+  AliHelix mhelix(x2,p2,sign);
+  
+  //
+  //
+  //
+  //find intersection linear
+  //
+  Double_t distance1, distance2;
+  Double_t phase[2][2],radius[2];
+  Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+  Double_t delta1=10000,delta2=10000;    
+  if (points>0){
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+  }
+  else{
+    fInvMass=-1;
+    return;
+  }
+  if (points==2){    
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+  }
+  distance1 = TMath::Min(delta1,delta2);
+  //
+  //find intersection parabolic
+  //
+  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+  delta1=10000,delta2=10000;  
+  
+  if (points>0){
+    dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+  }
+  if (points==2){    
+    dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+  }
+  
+  distance2 = TMath::Min(delta1,delta2);
+  //
+  if (delta1<delta2){
+    //get V0 info
+    dhelix1.Evaluate(phase[0][0],fMCXr);
+    dhelix1.GetMomentum(phase[0][0],fMCPdr);
+    mhelix.GetMomentum(phase[0][1],fMCPm);
+    dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
+    fMCRr = TMath::Sqrt(radius[0]);
+  }
+  else{
+    dhelix1.Evaluate(phase[1][0],fMCXr);
+    dhelix1.GetMomentum(phase[1][0],fMCPdr);
+    mhelix.GetMomentum(phase[1][1],fMCPm);
+    dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
+    fMCRr = TMath::Sqrt(radius[1]);
+  }     
+  //            
+  //
+  fMCDist1 = TMath::Sqrt(distance1);
+  fMCDist2 = TMath::Sqrt(distance2);      
+  //
+  //
+  // 
+  Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
+  //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
+  Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
+  Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
+  Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
+  vnorm2 = TMath::Sqrt(vnorm2);
+  Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
+  Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
+  pnorm2 = TMath::Sqrt(pnorm2);
+  //
+  fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
+  fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
+  fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
+  Float_t mass1 = fMCd.fMass;
+  Float_t mass2 = fMCm.fMass;
+
+  
+  Float_t e1    = TMath::Sqrt(mass1*mass1+
+                             fMCPdr[0]*fMCPdr[0]+
+                             fMCPdr[1]*fMCPdr[1]+
+                             fMCPdr[2]*fMCPdr[2]);
+  Float_t e2    = TMath::Sqrt(mass2*mass2+
+                             fMCPm[0]*fMCPm[0]+
+                             fMCPm[1]*fMCPm[1]+
+                             fMCPm[2]*fMCPm[2]);
+  
+  fInvMass =  
+    (fMCPm[0]+fMCPdr[0])*(fMCPm[0]+fMCPdr[0])+
+    (fMCPm[1]+fMCPdr[1])*(fMCPm[1]+fMCPdr[1])+
+    (fMCPm[2]+fMCPdr[2])*(fMCPm[2]+fMCPdr[2]);
+  
+  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
+
+    
+}
+
+
+
+
+
+/////////////////////////////////////////////////////////////////////////////////
+void AliGenKinkInfo::Update()
+{
+  fMCPd[0] = fMCd.fParticle.Px();
+  fMCPd[1] = fMCd.fParticle.Py();
+  fMCPd[2] = fMCd.fParticle.Pz();
+  fMCPd[3] = fMCd.fParticle.P();
+  //
+  fMCX[0]  = fMCd.fParticle.Vx();
+  fMCX[1]  = fMCd.fParticle.Vy();
+  fMCX[2]  = fMCd.fParticle.Vz();
+  fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
+  //
+  fPdg[0]    = fMCd.fParticle.GetPdgCode();
+  fPdg[1]    = fMCm.fParticle.GetPdgCode();
+  //
+  fLab[0]    = fMCd.fParticle.GetUniqueID();
+  fLab[1]    = fMCm.fParticle.GetUniqueID();
+  //
+  //
+  //
+  Double_t x1[3],p1[3];
+  TParticle & pdaughter = fMCd.fParticle;
+  x1[0] = pdaughter.Vx();      
+  x1[1] = pdaughter.Vy();
+  x1[2] = pdaughter.Vz();
+  p1[0] = pdaughter.Px();
+  p1[1] = pdaughter.Py();
+  p1[2] = pdaughter.Pz();
+  Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
+  AliHelix dhelix1(x1,p1,sign);
+  //
+  //
+  Double_t x2[3],p2[3];            
+  //
+  TParticle & pmother = fMCm.fParticle;
+  x2[0] = pmother.Vx();      
+  x2[1] = pmother.Vy();      
+  x2[2] = pmother.Vz();      
+  p2[0] = pmother.Px();
+  p2[1] = pmother.Py();
+  p2[2] = pmother.Pz();
+  //
+  AliTrackReference & pdecay = fMCm.fTRdecay;
+  x2[0] = pdecay.X();      
+  x2[1] = pdecay.Y();      
+  x2[2] = pdecay.Z();      
+  p2[0] = pdecay.Px();
+  p2[1] = pdecay.Py();
+  p2[2] = pdecay.Pz();
+  //
+  sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
+  AliHelix mhelix(x2,p2,sign);
+  
+  //
+  //
+  //
+  //find intersection linear
+  //
+  Double_t distance1, distance2;
+  Double_t phase[2][2],radius[2];
+  Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+  Double_t delta1=10000,delta2=10000;    
+  if (points>0){
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+  }
+  if (points==2){    
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+  }
+  distance1 = TMath::Min(delta1,delta2);
+  //
+  //find intersection parabolic
+  //
+  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
+  delta1=10000,delta2=10000;  
+  
+  if (points>0){
+    dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
+  }
+  if (points==2){    
+    dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
+  }
+  
+  distance2 = TMath::Min(delta1,delta2);
+  //
+  if (delta1<delta2){
+    //get V0 info
+    dhelix1.Evaluate(phase[0][0],fMCXr);
+    dhelix1.GetMomentum(phase[0][0],fMCPdr);
+    mhelix.GetMomentum(phase[0][1],fMCPm);
+    dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
+    fMCRr = TMath::Sqrt(radius[0]);
+  }
+  else{
+    dhelix1.Evaluate(phase[1][0],fMCXr);
+    dhelix1.GetMomentum(phase[1][0],fMCPdr);
+    mhelix.GetMomentum(phase[1][1],fMCPm);
+    dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
+    fMCRr = TMath::Sqrt(radius[1]);
+  }     
+  //            
+  //
+  fMCDist1 = TMath::Sqrt(distance1);
+  fMCDist2 = TMath::Sqrt(distance2);      
+    
 }
 
 
+Float_t AliGenKinkInfo::GetQt(){
+  //
+  //
+  Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
+  return TMath::Sin(fMCAngle[2])*momentumd;
+}
+
 
 
   
@@ -439,6 +807,7 @@ Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
   fStack = fLoader->Stack();
   //
   fLoader->LoadTrackRefs();
+  fLoader->LoadHits();
   fTreeTR = fLoader->TreeTR();
   //
   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
@@ -503,6 +872,12 @@ Int_t AliGenInfoMaker::Exec()
     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
     if (TreeKLoop()>0) return 1;
     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
+    //
+    if (BuildKinkInfo()>0) return 1;
+    if (BuildV0Info()>0) return 1;
+    //if (BuildHitLines()>0) return 1;
+    if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
+    //
     for (UInt_t i = 0; i<fNParticles; i++) {
       if (fGenInfo[i]) delete fGenInfo[i]; 
     }
@@ -529,10 +904,34 @@ void AliGenInfoMaker::CreateTreeGenTracks()
   }
   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");  
   AliMCInfo * info = new AliMCInfo;
-  //
   fTreeGenTracks->Branch("MC","AliMCInfo",&info);
   delete info; 
+  //
+  AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
+  fTreeKinks = new TTree("genKinksTree","genKinksTree");  
+  fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
+  delete kinkinfo;
+  //
+  AliGenV0Info *v0info = new AliGenV0Info;
+  fTreeV0 = new TTree("genV0Tree","genV0Tree");  
+  fTreeV0->Branch("MC","AliGenV0Info",&v0info);
+  delete v0info;
+  //
+  //
+  AliPointsMI * points0 = new AliPointsMI;  
+  AliPointsMI * points1 = new AliPointsMI;  
+  AliPointsMI * points2 = new AliPointsMI;  
+  fTreeHitLines = new TTree("HitLines","HitLines");
+  fTreeHitLines->Branch("TPC.","AliPointsMI",&points0,32000,10);
+  fTreeHitLines->Branch("TRD.","AliPointsMI",&points1,32000,10);
+  fTreeHitLines->Branch("ITS.","AliPointsMI",&points2,32000,10);
+  Int_t label=0;
+  fTreeHitLines->Branch("Label",&label,"label/I");
+  //
   fTreeGenTracks->AutoSave();
+  fTreeKinks->AutoSave();
+  fTreeV0->AutoSave();
+  fTreeHitLines->AutoSave();
 }
 ////////////////////////////////////////////////////////////////////////
 void AliGenInfoMaker::CloseOutputFile() 
@@ -544,6 +943,11 @@ void AliGenInfoMaker::CloseOutputFile()
   fFileGenTracks->cd();
   fTreeGenTracks->Write();  
   delete fTreeGenTracks;
+  fTreeKinks->Write();
+  delete fTreeKinks;
+  fTreeV0->Write();
+  delete fTreeV0;
+
   fFileGenTracks->Close();
   delete fFileGenTracks;
   return;
@@ -597,7 +1001,7 @@ Int_t AliGenInfoMaker::TreeKLoop()
     info->Update();
     //////////////////////////////////////////////////////////////////////    
     br->SetAddress(&info);    
-    fTreeGenTracks->Fill();
+    fTreeGenTracks->Fill();    
   }
   fTreeGenTracks->AutoSave();
   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
@@ -605,6 +1009,218 @@ Int_t AliGenInfoMaker::TreeKLoop()
 }
 
 
+////////////////////////////////////////////////////////////////////////////////////
+Int_t  AliGenInfoMaker::BuildKinkInfo()
+{
+  //
+  // Build tree with Kink Information
+  //
+  AliStack * stack = fStack;
+  if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
+  
+  if (fDebug > 0) {
+    cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
+       <<fEventNr<<endl;
+  }  
+  //  Int_t  ipdg = 0;
+  //TParticlePDG *ppdg = 0;
+  // not all generators give primary vertex position. Take the vertex
+  // of the particle 0 as primary vertex.
+  TDatabasePDG  pdg; //get pdg table  
+  //thank you very much root for this
+  TBranch * br = fTreeKinks->GetBranch("MC");
+  //
+  AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
+  //
+  for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+    // load only particles with TR
+    AliMCInfo * info = GetInfo(iParticle);
+    if (!info) continue;
+    if (info->fCharge==0) continue;  
+    if (info->fTRdecay.P()<0.13) continue;  //momenta cut 
+    if (info->fTRdecay.R()>500)  continue;  //R cut - decay outside barrel
+    TParticle & particle = info->fParticle;
+    Int_t first = particle.GetDaughter(0);
+    if (first ==0) continue;
+
+    Int_t last = particle.GetDaughter(1);
+    if (last ==0) last=first;
+    AliMCInfo * info2 =  0;
+    AliMCInfo * dinfo =  0;
+    
+    for (Int_t id2=first;id2<=last;id2++){
+      info2 = GetInfo(id2);
+      if (!info2) continue;
+      TParticle &p2 = info2->fParticle;
+      Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
+      if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
+      if (!(p2.GetPDG())) continue;
+      dinfo =info2;
+     
+      kinkinfo->fMCm = (*info);
+      kinkinfo->fMCd = (*dinfo);
+      if (kinkinfo->fMCm.fParticle.GetPDG()==0 ||  kinkinfo->fMCd.fParticle.GetPDG()==0) continue;
+      kinkinfo->Update();
+      br->SetAddress(&kinkinfo);    
+      fTreeKinks->Fill();
+    }
+    /*
+    if (dinfo){
+      kinkinfo->fMCm = (*info);
+      kinkinfo->fMCd = (*dinfo);
+      kinkinfo->Update();
+      br->SetAddress(&kinkinfo);    
+      fTreeKinks->Fill();
+    }
+    */
+  }
+  fTreeGenTracks->AutoSave();
+  if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
+  return 0;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////
+Int_t  AliGenInfoMaker::BuildV0Info()
+{
+  //
+  // Build tree with V0 Information
+  //
+  AliStack * stack = fStack;
+  if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
+  
+  if (fDebug > 0) {
+    cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
+       <<fEventNr<<endl;
+  }  
+  //  Int_t  ipdg = 0;
+  //TParticlePDG *ppdg = 0;
+  // not all generators give primary vertex position. Take the vertex
+  // of the particle 0 as primary vertex.
+  TDatabasePDG  pdg; //get pdg table  
+  //thank you very much root for this
+  TBranch * br = fTreeV0->GetBranch("MC");
+  //
+  AliGenV0Info * v0info = new AliGenV0Info;
+  //
+  //
+  for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+    // load only particles with TR
+    AliMCInfo * info = GetInfo(iParticle);
+    if (!info) continue;
+    if (info->fCharge==0) continue;  
+    //
+    //
+    TParticle & particle = info->fParticle;
+    Int_t mother = particle.GetMother(0);
+    if (mother <=0) continue;
+    //
+    TParticle * motherparticle = stack->Particle(mother);
+    if (!motherparticle) continue;
+    //
+    Int_t last = motherparticle->GetDaughter(1);
+    if (last==(int)iParticle) continue;
+    AliMCInfo * info2 =  info;
+    AliMCInfo * dinfo =  GetInfo(last);
+    if (!dinfo) continue;
+    if (!dinfo->fParticle.GetPDG()) continue;
+    if (!info2->fParticle.GetPDG()) continue;
+    if (dinfo){
+      v0info->fMCm = (*info);
+      v0info->fMCd = (*dinfo);
+      v0info->Update(fVPrim);
+      br->SetAddress(&v0info);    
+      fTreeV0->Fill();
+    }
+  }
+  fTreeV0->AutoSave();
+  if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
+  return 0;
+}
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::BuildHitLines()
+{
+
+//
+// open the file with treeK
+// loop over all entries there and save information about some tracks
+//
+
+  AliStack * stack = fStack;
+  if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
+  
+  if (fDebug > 0) {
+    cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
+       <<fEventNr<<endl;
+  }  
+  Int_t  ipdg = 0;
+  // TParticlePDG *ppdg = 0;
+  // not all generators give primary vertex position. Take the vertex
+  // of the particle 0 as primary vertex.
+  TDatabasePDG  pdg; //get pdg table  
+  //thank you very much root for this
+  AliPointsMI *tpcp = new AliPointsMI;
+  AliPointsMI *trdp = new AliPointsMI;
+  AliPointsMI *itsp = new AliPointsMI;
+  Int_t label =0;
+  //
+  TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
+  TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");  
+  TBranch * brits = fTreeHitLines->GetBranch("ITS.");
+  TBranch * brlabel = fTreeHitLines->GetBranch("Label");
+  brlabel->SetAddress(&label);
+  brtpc->SetAddress(&tpcp);
+  brtrd->SetAddress(&trdp);
+  brits->SetAddress(&itsp);
+  //
+  AliDetector *dtpc = gAlice->GetDetector("TPC");
+  AliDetector *dtrd = gAlice->GetDetector("TRD");
+  AliDetector *dits = gAlice->GetDetector("ITS");
+  for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
+    // load only particles with TR
+    AliMCInfo * info = GetInfo(iParticle);
+    if (!info) continue;
+    Int_t primpart = info->fPrimPart;
+    ipdg = info->fParticle.GetPdgCode();
+    label = iParticle;
+    //
+    gAlice->ResetHits();
+    tpcp->Reset();
+    itsp->Reset();
+    trdp->Reset();
+    tpcp->fLabel1 = ipdg;
+    trdp->fLabel1 = ipdg;
+    itsp->fLabel1 = ipdg;
+    if (dtpc->TreeH()->GetEvent(primpart)){
+      dtpc->LoadPoints(primpart);
+      tpcp->Reset(dtpc,iParticle);
+    }
+    if (dtrd->TreeH()->GetEvent(primpart)){
+      dtrd->LoadPoints(primpart);
+      trdp->Reset(dtrd,iParticle);
+    }
+    if (dits->TreeH()->GetEvent(primpart)){
+      dits->LoadPoints(primpart);
+      itsp->Reset(dits,iParticle);
+    }    
+    //    
+    fTreeHitLines->Fill();
+    dtpc->ResetPoints();
+    dtrd->ResetPoints();
+    dits->ResetPoints();
+  }
+  delete tpcp;
+  delete trdp;
+  delete itsp;
+  fTreeHitLines->AutoSave();
+  if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
+  return 0;
+}
 
 
 ////////////////////////////////////////////////////////////////////////
@@ -684,11 +1300,13 @@ Int_t AliGenInfoMaker::TreeTRLoop()
   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
   TClonesArray* ITSArrayTR = new TClonesArray("AliTrackReference");
   TClonesArray* TRDArrayTR = new TClonesArray("AliTrackReference");
+  TClonesArray* TOFArrayTR = new TClonesArray("AliTrackReference");
   TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
   //
   if (treeTR->GetBranch("TPC"))    treeTR->GetBranch("TPC")->SetAddress(&TPCArrayTR);
   if (treeTR->GetBranch("ITS"))    treeTR->GetBranch("ITS")->SetAddress(&ITSArrayTR);
   if (treeTR->GetBranch("TRD"))    treeTR->GetBranch("TRD")->SetAddress(&TRDArrayTR);
+  if (treeTR->GetBranch("TOF"))    treeTR->GetBranch("TOF")->SetAddress(&TOFArrayTR);
   if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR);
   //
   //
@@ -701,11 +1319,12 @@ Int_t AliGenInfoMaker::TreeTRLoop()
     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);            
       //
-      if (trackRef->Pt()<fgTPCPtCut) continue;
+      if (trackRef->P()<fgTPCPtCut) continue;
       Int_t label = trackRef->GetTrack();      
       AliMCInfo * info = GetInfo(label);
       if (!info) info = MakeInfo(label);
       if (!info) continue;
+      info->fPrimPart =  iPrimPart;
       TClonesArray & arr = *(info->fTPCReferences);
       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
     }
@@ -715,12 +1334,13 @@ Int_t AliGenInfoMaker::TreeTRLoop()
     for (Int_t iTrackRef = 0; iTrackRef < ITSArrayTR->GetEntriesFast(); iTrackRef++) {
       AliTrackReference *trackRef = (AliTrackReference*)ITSArrayTR->At(iTrackRef);            
       //
-      if (trackRef->Pt()<fgTPCPtCut) continue;
+      if (trackRef->P()<fgTPCPtCut) continue;
       Int_t label = trackRef->GetTrack();      
       AliMCInfo * info = GetInfo(label);
       if ( (!info) && trackRef->Pt()<fgITSPtCut) continue;
       if (!info) info = MakeInfo(label);
       if (!info) continue;
+      info->fPrimPart =  iPrimPart;
       TClonesArray & arr = *(info->fITSReferences);
       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
     }
@@ -730,16 +1350,34 @@ Int_t AliGenInfoMaker::TreeTRLoop()
     for (Int_t iTrackRef = 0; iTrackRef < TRDArrayTR->GetEntriesFast(); iTrackRef++) {
       AliTrackReference *trackRef = (AliTrackReference*)TRDArrayTR->At(iTrackRef);            
       //
-      if (trackRef->Pt()<fgTPCPtCut) continue;
+      if (trackRef->P()<fgTPCPtCut) continue;
       Int_t label = trackRef->GetTrack();      
       AliMCInfo * info = GetInfo(label);
       if ( (!info) && trackRef->Pt()<fgTRDPtCut) continue;
       if (!info) info = MakeInfo(label);
       if (!info) continue;
+      info->fPrimPart =  iPrimPart;
       TClonesArray & arr = *(info->fTRDReferences);
       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
     }
     //
+    // Loop over TOF references
+    //
+    for (Int_t iTrackRef = 0; iTrackRef < TOFArrayTR->GetEntriesFast(); iTrackRef++) {
+      AliTrackReference *trackRef = (AliTrackReference*)TOFArrayTR->At(iTrackRef);            
+      Int_t label = trackRef->GetTrack();      
+      AliMCInfo * info = GetInfo(label);
+      if (!info){
+       if (trackRef->Pt()<fgTPCPtCut) continue;
+       if ( (!info) && trackRef->Pt()<fgTOFPtCut) continue;
+      }
+      if (!info) info = MakeInfo(label);
+      if (!info) continue;
+      info->fPrimPart =  iPrimPart;
+      TClonesArray & arr = *(info->fTOFReferences);
+      new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
+    }
+    //
     // get dacay position
     for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
       AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);      
@@ -757,6 +1395,9 @@ Int_t AliGenInfoMaker::TreeTRLoop()
   delete  TPCArrayTR;
   TRDArrayTR->Delete();
   delete  TRDArrayTR;
+  TOFArrayTR->Delete();
+  delete  TOFArrayTR;
+
   ITSArrayTR->Delete();
   delete  ITSArrayTR;
   RunArrayTR->Delete();
@@ -1041,4 +1682,81 @@ TH1F* AliComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t draw
 }
 
 
+void   AliComparisonDraw::DrawFriend2D(const char * chx, const char *chy, const char* selection, TTree * tfriend)
+{
+  
+}
 
+
+void   AliComparisonDraw::GetPoints3D(const char * label, const char * chpoints, const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
+  //
+  //
+  if (!fPoints) fPoints = new TObjArray;
+  if (tpoints->GetIndex()==0) tpoints->BuildIndex("Label","Label");
+  AliPointsMI * points = new AliPointsMI;
+  TBranch * br = tpoints->GetBranch(chpoints);
+  br->SetAddress(&points);
+  Int_t npoints = fTree->Draw(label,selection);
+  Float_t xyz[30000];
+  rmin*=rmin;
+  for (Int_t i=0;i<npoints;i++){
+    Int_t index = (Int_t)fTree->GetV1()[i];
+    tpoints->GetEntryWithIndex(index,index);
+    if (points->fN<2) continue;
+    Int_t goodpoints=0;
+    for (Int_t i=0;i<points->fN; i++){
+      xyz[goodpoints*3]   = points->fX[i];
+      xyz[goodpoints*3+1] = points->fY[i];
+      xyz[goodpoints*3+2] = points->fZ[i];
+      if ( points->fX[i]*points->fX[i]+points->fY[i]*points->fY[i]>rmin) goodpoints++;
+    } 
+    TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz); 
+    //    marker->SeWidth(1); 
+    marker->SetMarkerColor(color);
+    marker->SetMarkerStyle(1);
+    fPoints->AddLast(marker);
+  }
+}
+
+void AliComparisonDraw::Draw3D(Int_t min, Int_t max){
+  if (!fPoints) return;
+  Int_t npoints = fPoints->GetEntriesFast();
+  max = TMath::Min(npoints,max);
+  min = TMath::Min(npoints,min);
+
+  for (Int_t i =min; i<max;i++){
+    TObject * obj = fPoints->At(i);
+    if (obj) obj->Draw();
+  }
+}
+
+void AliComparisonDraw:: SavePoints(const char* name)
+{
+  char fname[25];
+  sprintf(fname,"TH_%s.root",name);
+  TFile f(fname,"new");
+  fPoints->Write("Tracks",TObject::kSingleKey);
+  f.Close();
+  sprintf(fname,"TH%s.gif",name);
+  fCanvas->SaveAs(fname);
+}
+
+void   AliComparisonDraw::InitView(){
+  //
+  //
+  AliRunLoader* rl = AliRunLoader::Open();
+  rl->GetEvent(0); 
+  rl->CdGAFile(); 
+  //
+  fCanvas = new TCanvas("cdisplay", "Cluster display",0,0,700,730);
+  fView   = new TView(1);
+  fView->SetRange(-330,-360,-330, 360,360, 330);
+  fCanvas->Clear();
+  fCanvas->SetFillColor(1);
+  fCanvas->SetTheta(90.);
+  fCanvas->SetPhi(0.);
+  //
+  TGeometry *geom =(TGeometry*)gDirectory->Get("AliceGeom");
+  geom->Draw("same");
+
+}
index b80ee5805b86157484a14437814632665ee0644e..631f1385703344b185718c94ba17547cb1e74181 100644 (file)
@@ -45,6 +45,7 @@ public:
   AliTrackReference  fTrackRefOut;   // decay track reference saved in the output tree
   AliTrackReference  fTRdecay;       // track reference at decay point
   //
+  Int_t fPrimPart;                   // index of primary particle in TreeH
   TParticle fParticle;           // generated particle 
   Float_t   fMass;               // mass of the particle
   Float_t fCharge;               //
@@ -61,10 +62,13 @@ public:
   Float_t fPrim;               // theoretical dedx in tpc according particle momenta and mass
   digitRow fTPCRow;               // information about digits row pattern
   Int_t fNTPCRef;                 // tpc references counter
-  Int_t fNITSRef;                 // tpc references counter
+  Int_t fNITSRef;                 // ITS references counter
+  Int_t fNTRDRef;                  // TRD references counter
+  Int_t fNTOFRef;                  // TOF references counter
   TClonesArray * fTPCReferences;  //containner with all track references -in the TPC
   TClonesArray * fITSReferences;  //container with ITS references
   TClonesArray * fTRDReferences;  //container with TRD references  
+  TClonesArray * fTOFReferences;  //container with TRD references  
   //
   ClassDef(AliMCInfo,3)  // container for 
 };
@@ -74,9 +78,9 @@ ClassImp(AliMCInfo)
 
 class AliGenV0Info: public TObject {
 public:
-  AliMCInfo fMCd;      //info about daughter particle
-  AliMCInfo fMCm;      //info about mother particle
-  void Update();        // put some derived info to special field 
+  AliMCInfo fMCd;      //info about daughter particle - second particle for V0
+  AliMCInfo fMCm;      //info about mother particle   - first particle for V0
+  void Update(Float_t vertex[3]);        // put some derived info to special field 
   Double_t    fMCDist1;    //info about closest distance according closest MC - linear DCA
   Double_t    fMCDist2;    //info about closest distance parabolic DCA
   //
@@ -91,6 +95,12 @@ public:
   Double_t     fMCR;        //exact r position of the vertex
   Int_t        fPdg[2];   //pdg code of mother and daugter particles
   Int_t        fLab[2];   //MC label of the partecle
+  //
+  Double_t       fInvMass;  //reconstructed invariant mass -
+  Float_t        fPointAngleFi; //point angle fi
+  Float_t        fPointAngleTh; //point angle theta
+  Float_t        fPointAngle;   //point angle full
+  //
   ClassDef(AliGenV0Info,1)  // container for  
 };
 ClassImp(AliGenV0Info)
@@ -98,6 +108,33 @@ ClassImp(AliGenV0Info)
 
 
 
+class AliGenKinkInfo: public TObject {
+public:
+  AliMCInfo fMCd;      //info about daughter particle - second particle for V0
+  AliMCInfo fMCm;      //info about mother particle   - first particle for V0
+  void Update();        // put some derived info to special field 
+  Float_t GetQt();      //
+  Double_t    fMCDist1;    //info about closest distance according closest MC - linear DCA
+  Double_t    fMCDist2;    //info about closest distance parabolic DCA
+  //
+  Double_t     fMCPdr[3];    //momentum at vertex daughter  - according approx at DCA
+  Double_t     fMCPd[4];     //exact momentum from MC info
+  Double_t     fMCX[3];      //exact position of the vertex
+  Double_t     fMCXr[3];     //rec. position according helix
+  //
+  Double_t     fMCPm[3];    //momentum at the vertex mother
+  Double_t     fMCAngle[3]; //three angels
+  Double_t     fMCRr;       // rec position of the vertex 
+  Double_t     fMCR;        //exact r position of the vertex
+  Int_t        fPdg[2];   //pdg code of mother and daugter particles
+  Int_t        fLab[2];   //MC label of the partecle
+  ClassDef(AliGenKinkInfo,1)  // container for  
+};
+ClassImp(AliGenKinkInfo)
+
+
+
+
 ////////////////////////////////////////////////////////////////////////
 // 
 // Start of implementation of the class AliGenInfoMaker
@@ -119,6 +156,9 @@ public:
   Int_t TreeKLoop();
   Int_t TreeTRLoop();
   Int_t TreeDLoop();
+  Int_t BuildKinkInfo();  // build information about MC kinks
+  Int_t BuildV0Info();  // build information about MC kinks
+  Int_t BuildHitLines();  // build information about MC kinks
   //void  FillInfo(Int_t iParticle);
   void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
   void SetNEvents(Int_t i) {fNEvents = i;}
@@ -139,8 +179,11 @@ public:
   Int_t fNEvents;                 //! number of events to process
   Int_t fFirstEventNr;            //! first event to process
   UInt_t fNParticles;              //! number of particles in TreeK
-  TTree *fTreeGenTracks;          //! output tree with generated tracks
-  char  fFnRes[1000];                   //! output file name with stored tracks
+  TTree * fTreeGenTracks;          //! output tree with generated tracks
+  TTree * fTreeKinks;             //!  output tree with Kinks
+  TTree * fTreeV0;                //!  output tree with V0
+  TTree * fTreeHitLines;          //! tree with hit lines
+  char  fFnRes[1000];             //! output file name with stored tracks
   TFile *fFileGenTracks;          //! output file with stored fTreeGenTracks
   //
   AliRunLoader * fLoader;         //! pointer to the run loader
@@ -152,7 +195,7 @@ public:
   Int_t fNInfos;                  //! number of tracks with infos
   //
   AliTPCParam* fParamTPC;         //! AliTPCParam
-  Double_t fVPrim[3];             //! primary vertex position
+  Float_t fVPrim[3];             //! primary vertex position
                                   // the fVDist[3] contains size of the 3-vector
 
 private:
@@ -164,9 +207,10 @@ private:
   static const Int_t seedRow22 = 130;  // seedRow12 - shift
   static const Double_t kRaddeg = 180./3.14159265358979312;
   // 
-  static const Double_t fgTPCPtCut = 0.1; // do not store particles with generated pT less than this
+  static const Double_t fgTPCPtCut = 0.03; // do not store particles with generated pT less than this
   static const Double_t fgITSPtCut = 0.2; // do not store particles with generated pT less than this
   static const Double_t fgTRDPtCut = 0.2; // do not store particles with generated pT less than this
+  static const Double_t fgTOFPtCut = 0.15; // do not store particles with generated pT less than this
  
   ClassDef(AliGenInfoMaker,1)    // class which creates and fills tree with TPCGenTrack objects
 };
@@ -176,6 +220,8 @@ ClassImp(AliGenInfoMaker)
 
 class AliComparisonDraw: public TObject{
 public:
+  AliComparisonDraw(){fPoints=0; fView=0;}
+  void InitView();
   TH1F * DrawXY(const char * chx, const char *chy, const char* selection, 
                const char * quality,Int_t nbins, Float_t minx, Float_t maxx, 
                Float_t miny, Float_t maxy);
@@ -193,15 +239,40 @@ public:
   static TH1F*  CreateEffHisto(TH1F* hGen, TH1F* hRec);
   static TH1F*  CreateResHisto(TH2F* hRes2, TH1F **phMean, 
                                Bool_t drawBinFits = kTRUE,Bool_t overflowBinFits = kFALSE);
+  void   DrawFriend2D(const char * chx, const char *chy, const char* selection, TTree * tfriend);
+  void   GetPoints3D(const char * label, const char * chpoints, const char* selection, TTree * tpoints, Int_t color=6, Float_t rmin=4.);
+  void   Draw3D(Int_t min=0, Int_t max = 10000);
+  void SavePoints(const char* name);
  public: 
   TTree * fTree;
   TH1F  * fRes;  //temporary file
   TH1F  * fMean;  //temporary file
+  TView * fView;  //3D view
+  TCanvas *fCanvas; //canvas
+  TObjArray *fPoints;
   ClassDef(AliComparisonDraw,1)
 };
 ClassImp(AliComparisonDraw)
 
 
+class AliPointsMI: public TObject{
+ public:
+  AliPointsMI();
+  AliPointsMI(Int_t n, Float_t *x,Float_t *y, Float_t *z);
+  void Reset();
+  void Reset(AliDetector * det, Int_t particle);  //load points for given particle
+  ~AliPointsMI();
+  Int_t   fN;  //number of points;
+  Float_t *fX; //[fN] pointer to x
+  Float_t *fY; //[fN] pointer to y
+  Float_t *fZ; //[fN] pointer to Z  
+  Int_t   fCapacity; //!allocated size of the x,y,x
+  Int_t   fLabel0; //label
+  Int_t   fLabel1; //label
+  ClassDef(AliPointsMI,1)
+};
+ClassImp(AliPointsMI)
+
 
 AliTPCParam * GetTPCParam();
 Float_t TPCBetheBloch(Float_t bg);
similarity index 90%
rename from TPC/AliHelix.cxx
rename to STEER/AliHelix.cxx
index fd19f659586ba134556c050ca9767da38f5b0a79..78f250c26c0422a5c22e2101a1089d0c7df043b9 100644 (file)
@@ -23,6 +23,7 @@
 
 #include "AliHelix.h"
 #include "AliKalmanTrack.h"
+#include "AliExternalTrackParam.h"
 #include "TMath.h"
 ClassImp(AliHelix)
 
@@ -79,6 +80,44 @@ AliHelix::AliHelix(const AliKalmanTrack &t)
   //fHelix[0]+=TMath::Cos(fHelix[2])/fHelix[4];  
 }
 
+
+AliHelix::AliHelix(const AliExternalTrackParam &t)
+{
+  //
+  // 
+  Double_t alpha,x,cs,sn;
+  const Double_t *param =t.GetParameter(); 
+  for (Int_t i=0;i<5;i++) fHelix[i]=param[i]; 
+  x = t.X();
+  alpha=t.Alpha();
+  //
+  //circle parameters
+  fHelix[4]=fHelix[4]/AliKalmanTrack::GetConvConst();    // C
+  cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+
+  Double_t xc, yc, rc;
+  rc  =  1/fHelix[4];
+  xc  =  x-fHelix[2]*rc;
+  yc  =  fHelix[0]+TMath::Sqrt(1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4])/fHelix[4];
+  
+  fHelix[6] = xc*cs - yc*sn;
+  fHelix[7] = xc*sn + yc*cs;
+  fHelix[8] =  TMath::Abs(rc);
+  //
+  //
+  fHelix[5]=x*cs - fHelix[0]*sn;            // x0
+  fHelix[0]=x*sn + fHelix[0]*cs;            // y0
+  //fHelix[1]=                               // z0
+  fHelix[2]=TMath::ASin(fHelix[2]) + alpha; // phi0
+  //fHelix[3]=                               // tgl
+  //
+  //
+  fHelix[5]   = fHelix[6];
+  fHelix[0]   = fHelix[7];
+  //fHelix[5]-=TMath::Sin(fHelix[2])/fHelix[4]; 
+  //fHelix[0]+=TMath::Cos(fHelix[2])/fHelix[4];  
+}
+
 AliHelix::AliHelix(Double_t x[3], Double_t p[3], Double_t charge, Double_t conversion)
 {
   //
@@ -150,11 +189,19 @@ void   AliHelix::GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[
   Double_t norm2  = TMath::Sqrt(norm2r+g2[2]*g2[2]);
   norm2r         = TMath::Sqrt(norm2r);
   //
-  angle[0]  = TMath::ACos((g1[0]*g2[0]+g1[1]*g2[1])/(norm1r*norm2r));   // angle in phi projection
-  angle[1]  = TMath::ACos(((norm1r*norm2r)+g1[2]*g2[2])/(norm1*norm2)); // angle in rz  projection
-  angle[2]  = TMath::ACos((g1[0]*g2[0]+g1[1]*g2[1]+g1[2]*g2[2])/(norm1*norm2)); //3D angle
+  angle[0]  = (g1[0]*g2[0]+g1[1]*g2[1])/(norm1r*norm2r);   // angle in phi projection
+  if (TMath::Abs(angle[0])<1.) angle[0] = TMath::ACos(angle[0]);
+  else angle[0]=0;
+  //
+  angle[1]  = ((norm1r*norm2r)+g1[2]*g2[2])/(norm1*norm2); // angle in rz  projection
+  if (TMath::Abs(angle[1])<1.) angle[1] = TMath::ACos(angle[1]);
+  else angle[1]=0;
 
-    
+  angle[2]  = (g1[0]*g2[0]+g1[1]*g2[1]+g1[2]*g2[2])/(norm1*norm2); //3D angle
+  if (TMath::Abs(angle[2])<1.) angle[2] = TMath::ACos(angle[2]);
+  else angle[2]=0;
+
+  
   
 
 }
@@ -247,6 +294,7 @@ Int_t    AliHelix::GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Doubl
   Double_t  c2[3] = {h.fHelix[5]-fHelix[5],h.fHelix[0]-fHelix[0],h.fHelix[8]};
 
   Double_t d  = TMath::Sqrt(c2[0]*c2[0]+c2[1]*c2[1]); 
+  if (d<0.000000000001) return 0;
   //
   Double_t x0[2];
   Double_t y0[2];
@@ -255,11 +303,11 @@ Int_t    AliHelix::GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Doubl
     if (d>=(c1[2]+c2[2]+cut)) return 0;
     x0[0] = (d+c1[2]-c2[2])*c2[0]/(2*d)+ fHelix[5];
     y0[0] = (d+c1[2]-c2[2])*c2[1]/(2*d)+ fHelix[0];
-    return 0;
-//      phase[0][0] = GetPhase(x0[0],y0[0]);
-//      phase[0][1] = h.GetPhase(x0[0],y0[0]);
-//      ri[0] = x0[0]*x0[0]+y0[0]*y0[0];
-//      return 1;
+    //    return 0;
+    phase[0][0] = GetPhase(x0[0],y0[0]);
+    phase[0][1] = h.GetPhase(x0[0],y0[0]);
+    ri[0] = x0[0]*x0[0]+y0[0]*y0[0];
+    return 1;
   }
   if ( (d+c2[2])<c1[2]){
     if ( (d+c2[2])+cut<c1[2]) return 0;
similarity index 97%
rename from TPC/AliHelix.h
rename to STEER/AliHelix.h
index 35a0edd015f2ff78fb82767a4a04e46c6c835f32..de08f82e69985bf517217297be545d2085513382 100644 (file)
 
 class AliCluster;
 class AliKalmanTrack;
+class AliExternalTrackParam;
 
 class AliHelix : public TObject {
 public:
   AliHelix();
   AliHelix(const AliHelix &t);
   AliHelix(const AliKalmanTrack &t);
+  AliHelix(const AliExternalTrackParam &t);
   AliHelix(Double_t x[3], Double_t p[3], Double_t charge=1, Double_t conversion=0.);
   virtual ~AliHelix(){};
   inline void Evaluate(Double_t t, Double_t r[3]);
diff --git a/STEER/AliTrackParam.cxx b/STEER/AliTrackParam.cxx
new file mode 100644 (file)
index 0000000..cf9ae31
--- /dev/null
@@ -0,0 +1,240 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+// base class for track parameters                                           //
+//                                                                           //
+// The idea of having a separate class for the track parametrisation and     //
+// not including it in the track classes itself is to not duplicate code.    //
+// Track classes which use the same parametrisation can include a track      //
+// parameters object for their common parametrisation. The code, e.g. for    //
+// the propagation, does not need to be duplicated.                          //
+//                                                                           //
+// The AliTrackParam and its derived classes                                 //
+// - know the current track parameters and their covariance matrix as well   //
+//   as the local x coordinate and azimuthal angle alpha of the current      //
+//   parametrisation                                                         //
+// - can rotate their local coordinate system                                //
+// - can create a parametrisation in external format from itself             //
+// - can propagate through material or vacuum                                //
+// - can calculate a chi^2 for a cluster                                     //
+// - can update the parameters using the position and covariance of a        //
+//   cluster                                                                 //
+//                                                                           //
+// In addition some methods to get quantities useful for analysis, like      //
+// the momentum, are implemented.                                            //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+
+#include "AliTrackParam.h"
+#include "AliExternalTrackParam.h"
+#include "AliRunLoader.h"
+#include "AliRun.h"
+#include "AliMagF.h"
+
+
+Double_t FastMath::fgFastAsin[20000];
+
+FastMath::FastMath()
+{
+  for (Int_t i=0;i<10000;i++){
+    fgFastAsin[2*i] = TMath::ASin(i/10000.);
+    fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
+  }
+}
+
+Double_t FastMath::FastAsin(Double_t x)
+{
+  if (x>0){
+    Int_t index = int(x*10000);
+    return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
+  }
+  x*=-1;
+  Int_t index = int(x*10000);
+  return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
+}
+
+FastMath gFastMath;
+
+
+
+ClassImp(AliTrackParam)
+
+
+
+
+//_____________________________________________________________________________
+Bool_t AliTrackParam::RotateAndPropagateTo(Double_t alpha, Double_t x, 
+                                          Double_t* length)
+{
+// Rotate the reference axis for the parametrisation to the given angle and
+// propagate the track parameters to the given x coordinate assuming vacuum.
+// If length is not NULL, the change of track length is added to it.
+
+  if (!RotateTo(alpha)) return kFALSE;
+  if (!PropagateTo(x, length)) return kFALSE;
+  return kTRUE;
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::GetDsdx() const
+{
+// get the change of track length s per step in x: ds/dx
+
+  TVector3 x(TMath::Cos(Alpha()), TMath::Sin(Alpha()), 0);
+  TVector3 p = Momentum();
+  Double_t xp = x*p;
+  if (xp == 0) return 1.E6; 
+  return p.Mag() / xp;
+}
+
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::Phi() const
+{
+// get the azimuthal angre
+
+  return Momentum().Phi();
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::SigmaPhi() const
+{
+// get the error of the azimuthal angle
+
+  AliExternalTrackParam* param = CreateExternalParam();
+  Double_t result = param->SigmaPhi();
+  delete param;
+  return result;
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::Theta() const
+{
+// the the polar angle
+
+  return Momentum().Theta();
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::SigmaTheta() const
+{
+// get the error of the polar angle
+
+  AliExternalTrackParam* param = CreateExternalParam();
+  Double_t result = param->SigmaTheta();
+  delete param;
+  return result;
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::Eta() const
+{
+// get the pseudorapidity
+
+  return Momentum().Eta();
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::Px() const
+{
+// get the x component of the momentum
+
+  return Momentum().Px();
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::Py() const
+{
+// get the y component of the momentum
+
+  return Momentum().Py();
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::Pz() const
+{
+// get the z component of the momentum
+
+  return Momentum().Pz();
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::Pt() const
+{
+// get the transversal component of the momentum
+
+  return Momentum().Pt();
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::SigmaPt() const
+{
+// get the error of the transversal component of the momentum
+
+  AliExternalTrackParam* param = CreateExternalParam();
+  Double_t result = param->SigmaPt();
+  delete param;
+  return result;
+}
+
+//_____________________________________________________________________________
+Double_t AliTrackParam::P() const
+{
+// get the absolute momentum
+
+  return Momentum().Mag();
+}
+
+//_____________________________________________________________________________
+TVector3 AliTrackParam::Momentum() const
+{
+// get the momentum vector
+
+  AliExternalTrackParam* param = CreateExternalParam();
+  TVector3 result = param->Momentum();
+  delete param;
+  return result;
+}
+
+//_____________________________________________________________________________
+TVector3 AliTrackParam::Position() const
+{
+// get the current spatial position in global coordinates
+
+  Double_t sinAlpha = TMath::Sin(Alpha());
+  Double_t cosAlpha = TMath::Cos(Alpha());
+  return TVector3(X()*cosAlpha - Y()*sinAlpha, 
+                 X()*sinAlpha + Y()*cosAlpha, 
+                 Z());
+}
+
+//_____________________________________________________________________________
+TVector3 AliTrackParam::PositionAt(Double_t x) const
+{
+// get the spatial position at x in global coordinates
+
+  Double_t y;
+  Double_t z;
+  if (!GetProlongationAt(x, y, z)) return TVector3(0,0,0);
+  Double_t sinAlpha = TMath::Sin(Alpha());
+  Double_t cosAlpha = TMath::Cos(Alpha());
+  return TVector3(x*cosAlpha - y*sinAlpha, x*sinAlpha + y*cosAlpha, z);
+}
+
diff --git a/STEER/AliTrackParam.h b/STEER/AliTrackParam.h
new file mode 100644 (file)
index 0000000..561a939
--- /dev/null
@@ -0,0 +1,72 @@
+#ifndef ALITRACKPARAM_H
+#define ALITRACKPARAM_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+
+#include <TObject.h>
+#include <TVector3.h>
+
+class AliCluster;
+class AliExternalTrackParam;
+
+
+class FastMath {
+public:
+  FastMath();  
+  static Double_t FastAsin(Double_t x);   
+ private: 
+  static Double_t fgFastAsin[20000];
+};
+
+
+class AliTrackParam: public TObject {
+ public:
+  virtual const Double_t* GetParameter() const = 0;
+  virtual const Double_t* GetCovariance() const = 0;
+  virtual Double_t     X() const = 0;
+  virtual Double_t     Alpha() const = 0;
+  virtual AliExternalTrackParam* CreateExternalParam() const = 0;
+  virtual void         ResetCovariance(Double_t factor = 10.,
+                                      Bool_t clearOffDiagonal = kTRUE) = 0;
+
+  virtual Double_t     Y() const = 0;
+  virtual Double_t     Z() const = 0;
+
+  virtual Bool_t       PropagateTo(Double_t x, Double_t* length = NULL) = 0;
+  virtual Bool_t       RotateTo(Double_t alpha) = 0;
+  Bool_t               RotateBy(Double_t dAlpha)
+    {return RotateTo(Alpha() + dAlpha);};
+  virtual Bool_t       RotateAndPropagateTo(Double_t alpha, Double_t x, 
+                                           Double_t* length);
+  virtual Bool_t       CorrectForMaterial(Double_t dAngle2, 
+                                         Double_t dPrel) = 0;
+  virtual Bool_t       GetProlongationAt(Double_t x, Double_t& y, 
+                                        Double_t& z) const = 0;
+  virtual Double_t     GetXAtVertex(Double_t x = 0, Double_t y = 0) const = 0;
+  virtual Double_t     GetDsdx() const;
+
+  virtual Double_t     GetPredictedChi2(const AliCluster* cluster) = 0;
+  virtual Bool_t       Update(const AliCluster* cluster) = 0;
+
+  virtual Double_t     Phi() const;
+  virtual Double_t     SigmaPhi() const;
+  virtual Double_t     Theta() const;
+  virtual Double_t     SigmaTheta() const;
+  virtual Double_t     Eta() const;
+  virtual Double_t     Px() const;
+  virtual Double_t     Py() const;
+  virtual Double_t     Pz() const;
+  virtual Double_t     Pt() const;
+  virtual Double_t     SigmaPt() const;
+  virtual Double_t     P() const;
+  virtual TVector3     Momentum() const;
+  virtual TVector3     Position() const;
+  virtual TVector3     PositionAt(Double_t x) const;
+
+  ClassDef(AliTrackParam, 1)
+};
+
+#endif
index 8d03f879de59bbd0d1529efd505131f4b190e884..41423379ef41845b66e5054453ad5a85b317bbc1 100644 (file)
 #pragma link C++ class  AliESDcascade+;
 #pragma link C++ class  AliESDVertex+;
 #pragma link C++ class  AliESDpid+;
+#pragma link C++ class  AliESDkink+;
+#pragma link C++ class  AliESDV0MI+;
 
 #pragma link C++ class  AliKalmanTrack+;
+#pragma link C++ class  AliHelix+;
+#pragma link C++ class  AliExternalTrackParam+;
+#pragma link C++ class  AliTrackParam+;
 #pragma link C++ class  AliLog+;
 #endif
 
index a40b50132d5f3349686407e9cdfb26197105b8af..e1cc459abe3959cba114093bea22d33f5b0e07a1 100644 (file)
@@ -75,8 +75,3 @@
 #pragma link C++ class  AliVertexGenFile+;
 #pragma link C++ class  AliVertexer+;
 #endif
-
-
-
-
-
index 34c6ef409a7403ab9d053f3854ecb4343781e4ae..ecf1fd3f2fff664f5b7d84aa0a7c1e6a4fab306e 100644 (file)
@@ -2,9 +2,9 @@ SRCS = AliESD.cxx \
        AliESDtrack.cxx \
        AliESDMuonTrack.cxx AliESDPmdTrack.cxx AliESDHLTtrack.cxx \
        AliESDv0.cxx AliESDcascade.cxx AliESDVertex.cxx \
-       AliESDpid.cxx \
-       AliKalmanTrack.cxx \
-       AliLog.cxx
+       AliESDpid.cxx AliESDkink.cxx AliESDV0MI.cxx \
+       AliKalmanTrack.cxx AliHelix.cxx AliExternalTrackParam.cxx \
+       AliTrackParam.cxx AliLog.cxx
 
 HDRS:= $(SRCS:.cxx=.h) 
 
index 8e50623f63d394714b3cdc17ce08707ed40078e1..b8f6bd08139d19fe3b214b2995fad4b784f18b7a 100644 (file)
@@ -39,9 +39,13 @@ AliTPCtrack::AliTPCtrack(): AliKalmanTrack()
   fX = fP0 = fP1 = fP2 = fP3 = fP3 = fP4 = 0.0;
   fAlpha = fdEdx = 0.0;
   fNWrong = fNRotation = fNumber = 0;  // [SR, 01.04.2003]
+  for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
 }
 
 //_________________________________________________________________________
+
+
+
 AliTPCtrack::AliTPCtrack(UInt_t index, const Double_t xx[5],
 const Double_t cc[15], Double_t xref, Double_t alpha) : AliKalmanTrack() {
   //-----------------------------------------------------------------
@@ -73,6 +77,7 @@ const Double_t cc[15], Double_t xref, Double_t alpha) : AliKalmanTrack() {
   fRemoval    = 0;
   fTrackType  = 0;
   fLab2       = 0;
+  for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
 }
 
 //_____________________________________________________________________________
@@ -121,6 +126,7 @@ AliKalmanTrack(t) {
   fRemoval    = 0;
   fTrackType  = 0;
   fLab2       = 0;
+  for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
 }
 
 //_____________________________________________________________________________
@@ -131,6 +137,7 @@ AliTPCtrack::AliTPCtrack(const AliESDtrack& t) : AliKalmanTrack() {
   SetNumberOfClusters(t.GetTPCclusters(fIndex));
   SetLabel(t.GetLabel());
   SetMass(t.GetMass());
+  for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.GetKinkIndex(i);
 
   fdEdx  = t.GetTPCsignal();
   fAlpha = t.GetAlpha();
@@ -139,6 +146,8 @@ AliTPCtrack::AliTPCtrack(const AliESDtrack& t) : AliKalmanTrack() {
 
   //Conversion of the track parameters
   Double_t x,p[5]; t.GetExternalParameters(x,p);
+  Double_t c[15]; t.GetExternalCovariance(c);
+
   fX=x;    x=GetConvConst();
   fP0=p[0];
   fP1=p[1];
@@ -147,7 +156,6 @@ AliTPCtrack::AliTPCtrack(const AliESDtrack& t) : AliKalmanTrack() {
   fP2=fP4*fX - p[2];
 
   //Conversion of the covariance matrix
-  Double_t c[15]; t.GetExternalCovariance(c);
   c[10]/=x; c[11]/=x; c[12]/=x; c[13]/=x; c[14]/=x*x;
 
   Double_t c22=fX*fX*c[14] - 2*fX*c[12] + c[5];
@@ -206,7 +214,7 @@ AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : AliKalmanTrack(t) {
   fRemoval    = t.fRemoval ;
   fTrackType  = t.fTrackType;
   fLab2       = t.fLab2;
-
+  for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.fKinkIndexes[i];
 }
 //_____________________________________________________________________________
 
@@ -673,3 +681,65 @@ Double_t AliTPCtrack::GetD(Double_t x, Double_t y) const {
   delta -= TMath::Abs(r);
   return delta;  
 }
+
+//
+//
+
+void  AliTPCtrack::UpdatePoints()
+{
+  //--------------------------------------------------
+  //calculates first ,amx dens and last points
+  //--------------------------------------------------
+  Float_t density[160];
+  for (Int_t i=0;i<160;i++) density[i]=-1.;
+  fPoints[0]= 160;
+  fPoints[1] = -1;
+  //
+  Int_t ngood=0;
+  Int_t undeff=0;
+  Int_t nall =0;
+  Int_t range=20;
+  for (Int_t i=0;i<160;i++){
+    Int_t last = i-range;
+    if (nall<range) nall++;
+    if (last>=0){
+      if (fIndex[last]>0&& (fIndex[last]&0x8000)==0) ngood--;
+      if (fIndex[last]==-1) undeff--;
+    }
+    if (fIndex[i]>0&& (fIndex[i]&0x8000)==0)   ngood++;
+    if (fIndex[i]==-1) undeff++;
+    if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
+  }
+  Float_t maxdens=0;
+  Int_t indexmax =0;
+  for (Int_t i=0;i<160;i++){
+    if (density[i]<0) continue;
+    if (density[i]>maxdens){
+      maxdens=density[i];
+      indexmax=i;
+    }
+  }
+  //
+  //max dens point
+  fPoints[3] = maxdens;
+  fPoints[1] = indexmax;
+  //
+  // last point
+  for (Int_t i=indexmax;i<160;i++){
+    if (density[i]<0) continue;
+    if (density[i]<maxdens/2.) {
+      break;
+    }
+    fPoints[2]=i;
+  }
+  //
+  // first point
+  for (Int_t i=indexmax;i>0;i--){
+    if (density[i]<0) continue;
+    if (density[i]<maxdens/2.) {
+      break;
+    }
+    fPoints[0]=i;
+  }
+  //
+}
index adb8d6b2ce708ef28d6292e7817c91d40ee68257..199078ceea8a5c1b087023c067466a085788f864 100644 (file)
@@ -25,7 +25,7 @@
 #include <TMath.h>
 
 #include "AliTPCreco.h"
-
+#include "AliExternalTrackParam.h"
 class AliBarrelTrack;
 class AliESDtrack;
 
@@ -106,7 +106,8 @@ public:
   Int_t PropagateTo(Double_t xr,Double_t x0=28.94,Double_t rho=0.9e-3);
   Int_t Update(const AliCluster* c, Double_t chi2, UInt_t i);
   void ResetCovariance();
+  void UpdatePoints();//update points 
+  Float_t* GetPoints() {return fPoints;}
   //
   //  TClonesArray * fHelixIn;    //array of helixes  
   //TClonesArray * fHelixOut;   //array of helixes 
@@ -115,7 +116,10 @@ public:
   Float_t Density2(Int_t row0, Int_t row1); //calculate cluster density
   Double_t GetZat0() const;
   Double_t GetD(Double_t x=0, Double_t y=0) const;
-
+  AliExternalTrackParam & GetReference(){ return fReference;}
+  void UpdateReference(){ new (&fReference) AliExternalTrackParam(*this);}
+  Int_t  GetKinkIndex(Int_t i) const{ return fKinkIndexes[i];}
+  Int_t*  GetKinkIndexes() { return fKinkIndexes;}
 protected: 
   Double_t fX;              // X-coordinate of this track (reference plane)
   Double_t fAlpha;          // Rotation angle the local (TPC sector)
@@ -136,7 +140,7 @@ protected:
   Double_t fC40, fC41, fC42, fC43, fC44; // parameters
 
   Int_t fIndex[kMaxRow];       // indices of associated clusters 
-
+  Float_t fPoints[4];            //first, max dens row  end points of the track and max density
   //[SR, 01.04.2003]
   Int_t fNWrong;         // number of wrong clusters
   Int_t fNRotation;      // number of rotations
@@ -152,8 +156,9 @@ protected:
   Int_t fTrackType;       // track type - 0 - normal - 1 - kink -  2 -V0  3- double found
   Int_t fLab2;            // index of corresponding track (kink, V0, double)
   Int_t   fNShared;       // number of shared points 
+  AliExternalTrackParam   fReference; // track parameters at the middle of the chamber
   Float_t  fKinkPoint[12];      //radius, of kink,  dfi and dtheta
-
+  Int_t    fKinkIndexes[3];     // kink indexes - minus = mother + daughter
   ClassDef(AliTPCtrack,1)   // Time Projection Chamber reconstructed tracks
 };
 
index 2dda898c209b856ceb95b07bd3f8bbcea7afa465..1c4c60b77269ac3ce36608ea398131170c53b7a7 100644 (file)
@@ -51,6 +51,7 @@
 #include "TStopwatch.h"
 
 #include "AliTPCReconstructor.h"
+#include "AliESDkink.h"
 //
 
 ClassImp(AliTPCseed)
@@ -358,17 +359,35 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
     // write tracks to the event
     // store index of the track
     Int_t nseed=arr->GetEntriesFast();
+    //FindKinks(arr,fEvent);
     for (Int_t i=0; i<nseed; i++) {
       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
       if (!pt) continue; 
+      pt->UpdatePoints();
       pt->PropagateTo(fParam->GetInnerRadiusLow());
+      if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
+       AliESDtrack iotrack;
+       iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
+       iotrack.SetTPCPoints(pt->GetPoints());
+       iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+       //iotrack.SetTPCindex(i);
+       fEvent->AddTrack(&iotrack);
+       continue;
+      }
+       
       if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.55) {
        AliESDtrack iotrack;
-       iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);      
+       iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
+       iotrack.SetTPCPoints(pt->GetPoints());
        //iotrack.SetTPCindex(i);
+       iotrack.SetKinkIndexes(pt->GetKinkIndexes());
        fEvent->AddTrack(&iotrack);
        continue;
       } 
+      //
+      // short tracks  - maybe decays
+
       if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.70) {
        Int_t found,foundable,shared;
        pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
@@ -376,6 +395,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
          AliESDtrack iotrack;
          iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);    
          //iotrack.SetTPCindex(i);
+         iotrack.SetTPCPoints(pt->GetPoints());
+         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
          fEvent->AddTrack(&iotrack);
          continue;
        }
@@ -389,12 +410,46 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
        //
        AliESDtrack iotrack;
        iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);      
+       iotrack.SetTPCPoints(pt->GetPoints());
+       iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+       //iotrack.SetTPCindex(i);
+       fEvent->AddTrack(&iotrack);
+       continue;
+      }   
+      // short tracks  - secondaties
+      //
+      if ( (pt->GetNumberOfClusters()>30) ) {
+       Int_t found,foundable,shared;
+       pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
+       if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
+         AliESDtrack iotrack;
+         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);    
+         iotrack.SetTPCPoints(pt->GetPoints());
+         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+         //iotrack.SetTPCindex(i);
+         fEvent->AddTrack(&iotrack);
+         continue;
+       }
+      }       
+      
+      if ( (pt->GetNumberOfClusters()>15)) {
+       Int_t found,foundable,shared;
+       pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
+       if (found<15) continue;
+       if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
+       if (float(found)/float(foundable)<0.8) continue;
+       //
+       AliESDtrack iotrack;
+       iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);      
+       iotrack.SetTPCPoints(pt->GetPoints());
+       iotrack.SetKinkIndexes(pt->GetKinkIndexes());
        //iotrack.SetTPCindex(i);
        fEvent->AddTrack(&iotrack);
        continue;
       }   
     }
   }
+  printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
 }
 
 void AliTPCtrackerMI::WriteTracks(TTree * tree)
@@ -1842,7 +1897,31 @@ Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
     
   Int_t first = GetRowNumber(xt)-1;
-  for (Int_t nr= first; nr>=rf; nr-=step) {    
+  for (Int_t nr= first; nr>=rf; nr-=step) {
+    // update kink info
+    if (t.GetKinkIndexes()[0]>0){
+      for (Int_t i=0;i<3;i++){
+       Int_t index = t.GetKinkIndexes()[i];
+       if (index==0) break;
+       if (index<0) continue;
+       //
+       AliESDkink * kink = fEvent->GetKink(index-1);
+       if (!kink){
+         printf("PROBLEM\n");
+       }
+       else{
+         Int_t kinkrow = kink->fRow0+Int_t(1./(0.1+kink->fAngle[2]));
+         if (kinkrow==nr){
+           AliExternalTrackParam paramd(t);
+           kink->SetDaughter(paramd);
+           kink->Update();
+           kink->fStatus+=100;
+         }
+       }
+      }
+    }
+
+    if (nr==80) t.UpdateReference();
     if (nr<fInnerSec->GetNRows()) 
       fSectors = fInnerSec;
     else
@@ -1897,6 +1976,27 @@ Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
   if (first<0) first=0;
   for (Int_t nr=first; nr<=rf; nr++) {
     //if ( (t.GetSnp()<0.9))
+    if (t.GetKinkIndexes()[0]<0){
+      for (Int_t i=0;i<3;i++){
+       Int_t index = t.GetKinkIndexes()[i];
+       if (index==0) break;
+       if (index>0) continue;
+       index = TMath::Abs(index);
+       AliESDkink * kink = fEvent->GetKink(index-1);
+       if (!kink){
+         printf("PROBLEM\n");
+       }
+       else{
+         Int_t kinkrow = kink->fRow0-Int_t(1./(0.1+kink->fAngle[2]));
+         if (kinkrow==nr){
+           AliExternalTrackParam paramm(t);
+           kink->SetMother(paramm);
+           kink->Update();
+           kink->fStatus+=10;
+         }
+       }
+      }
+    }
     if (nr<fInnerSec->GetNRows()) 
       fSectors = fInnerSec;
     else
@@ -2226,6 +2326,88 @@ void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1,  Float_t fact
   }
 }
 
+
+void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t minimal)
+{
+
+  //Loop over all tracks and remove "overlaps"
+  //
+  //
+  UnsignClusters();
+  //
+  Int_t nseed = arr->GetEntriesFast();  
+  Float_t * quality = new Float_t[nseed];
+  Int_t   * indexes = new Int_t[nseed];
+  Int_t good =0;
+  //
+  //
+  for (Int_t i=0; i<nseed; i++) {
+    AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
+    if (!pt){
+      quality[i]=-1;
+      continue;
+    }
+    pt->UpdatePoints();    //select first last max dens points
+    Float_t * points = pt->GetPoints();
+    if (points[3]<0.8) quality[i] =-1;
+    //
+    quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
+  }
+  TMath::Sort(nseed,quality,indexes);
+  //
+  //
+  for (Int_t itrack=0; itrack<nseed; itrack++) {
+    Int_t trackindex = indexes[itrack];
+    AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);    
+    if (quality[trackindex]<0){
+      if (pt) {
+       delete arr->RemoveAt(trackindex);
+      }
+      else{
+       arr->RemoveAt(trackindex);
+      }
+      continue;
+    }
+    //
+    Int_t first = Int_t(pt->GetPoints()[0]);
+    Int_t last  = Int_t(pt->GetPoints()[2]);
+    Double_t factor = (pt->fBConstrain) ? factor1: factor2;
+    //
+    Int_t found,foundable,shared;
+    pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE);
+    Float_t sharedfactor = Float_t(shared)/Float_t(found);
+    //
+    if (Float_t(shared)/Float_t(found)>factor){
+      if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
+      delete arr->RemoveAt(trackindex);
+      continue;
+    }
+
+    if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){  //remove short tracks
+      if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
+      delete arr->RemoveAt(trackindex);
+      continue;
+    }
+
+    good++;
+    if (sharedfactor>0.4) continue;
+    for (Int_t i=first; i<last; i++) {
+      Int_t index=pt->GetClusterIndex2(i);
+      // if (index<0 || index&0x8000 ) continue;
+      if (index<0 || index&0x8000 ) continue;
+      AliTPCclusterMI *c= pt->fClusterPointer[i];        
+      if (!c) continue;
+      c->Use(10);  
+    }    
+  }
+  fNtracks = good;
+  if (fDebug>0){
+    Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
+  }
+  delete []quality;
+  delete []indexes;
+}
+
 void AliTPCtrackerMI::UnsignClusters() 
 {
   //
@@ -2490,17 +2672,20 @@ Int_t AliTPCtrackerMI::RefitInward(AliESD *event)
     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
     if (!seed) continue;
     seed->PropagateTo(fParam->GetInnerRadiusLow());
+    seed->UpdatePoints();
     AliESDtrack *esd=event->GetTrack(i);
     seed->CookdEdx(0.02,0.6);
     CookLabel(seed,0.1); //For comparison only
-    if (seed->GetNumberOfClusters()>60){
+    if (seed->GetNumberOfClusters()>15){
       esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit); 
+      esd->SetTPCPoints(seed->GetPoints());
       ntracks++;
     }
     else{
       //printf("problem\n");
     }
   }
+  //FindKinks(fSeeds,event);
   Info("RefitInward","Number of refitted tracks %d",ntracks);
   fEvent =0;
   //WriteTracks();
@@ -2523,14 +2708,18 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESD *event)
   for (Int_t i=0;i<nseed;i++){
     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
     if (!seed) continue;
+    if (seed->GetKinkIndex(0)<0)  UpdateKinkQualityM(seed);  // update quality informations for kinks
+    seed->UpdatePoints();
     AliESDtrack *esd=event->GetTrack(i);
     seed->CookdEdx(0.02,0.6);
     CookLabel(seed,0.1); //For comparison only
-    if (seed->GetNumberOfClusters()>60){
+    if (seed->GetNumberOfClusters()>15){
       esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
+      esd->SetTPCPoints(seed->GetPoints());
       ntracks++;
     }
   }
+  //FindKinks(fSeeds,event);
   Info("PropagateBack","Number of back propagated tracks %d",ntracks);
   fEvent =0;
   //WriteTracks();
@@ -2569,9 +2758,11 @@ void AliTPCtrackerMI::ReadSeeds(AliESD *event, Int_t direction)
   //  Int_t ntrk=0;  
   for (Int_t i=0; i<nentr; i++) {
     AliESDtrack *esd=event->GetTrack(i);
-    ULong_t status=esd->GetStatus();    
+    ULong_t status=esd->GetStatus();
+    if (!(status&AliESDtrack::kTPCin)) continue;
     AliTPCtrack t(*esd);
     AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
+    for (Int_t ikink=0;ikink<3;ikink++) seed->GetKinkIndexes()[ikink] = esd->GetKinkIndex(ikink);
     if ((status==AliESDtrack::kTPCin)&&(direction==1)) seed->ResetCovariance(); 
     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance();
     if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
@@ -2585,8 +2776,10 @@ void AliTPCtrackerMI::ReadSeeds(AliESD *event, Int_t direction)
       esd->GetExternalParameters(x,par1);
       Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
       Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
+      Double_t trdchi2=0;
+      if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
       //reset covariance if suspicious 
-      if ( (delta1>0.1) || (delta2>0.006))
+      if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
        seed->ResetCovariance();
     }
 
@@ -3723,6 +3916,338 @@ AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, F
   return seed;
 }
 
+void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
+{
+  //
+  //  find kinks
+  //
+  //
+  TObjArray *kinks= new TObjArray(10000);
+  Int_t nentries = array->GetEntriesFast();
+  AliHelix *helixes      = new AliHelix[nentries];
+  Int_t    *sign         = new Int_t[nentries];
+  Int_t    *nclusters    = new Int_t[nentries];
+  Float_t  *alpha        = new Float_t[nentries];
+  AliESDkink * kink      = new AliESDkink();
+  Int_t      * usage     = new Int_t[nentries];
+  Float_t  *zm             = new Float_t[nentries];
+  Float_t  *fim            = new Float_t[nentries];
+
+  //
+  //
+  for (Int_t i=0;i<nentries;i++){
+    sign[i]=0;
+    usage[i]=0;
+    AliTPCseed* track = (AliTPCseed*)array->At(i);    
+    if (!track) continue;
+    track->UpdatePoints();
+    if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
+      nclusters[i]=track->GetNumberOfClusters();
+      alpha[i] = track->GetAlpha();
+      new (&helixes[i]) AliHelix(*track);
+      sign[i] = (track->GetC()>0) ? -1:1;
+      Double_t x,y,z;
+      x=160;
+      if (track->GetProlongation(x,y,z)){
+       zm[i]  = z;
+       fim[i] = alpha[i]+TMath::ATan2(y,x);
+      }
+      else{
+       zm[i]  = track->GetZ();
+       fim[i] = alpha[i];
+      }
+    }
+  }
+  //
+  //
+  TStopwatch timer;
+  timer.Start();
+  Int_t ncandidates =0;
+  Int_t nall =0;
+  Int_t ntracks=0; 
+  Double_t phase[2][2],radius[2];
+  //
+  //
+  for (Int_t i =0;i<nentries;i++){
+    if (sign[i]==0) continue;
+    AliTPCseed * track0 = (AliTPCseed*)array->At(i);
+    ntracks++;
+    //
+    Double_t cradius0 = 40*40;
+    Double_t cradius1 = 270*270;
+    Double_t cdist1=8.;
+    Double_t cdist2=8.;
+    Double_t cdist3=0.55; 
+    for (Int_t j =i+1;j<nentries;j++){
+      nall++;
+      if (sign[j]*sign[i]<1) continue;
+      if ( (nclusters[i]+nclusters[j])>200) continue;
+      if ( (nclusters[i]+nclusters[j])<80) continue;
+      if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
+      if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
+      //AliTPCseed * track1 = (AliTPCseed*)array->At(j);  Double_t phase[2][2],radius[2];    
+      Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
+      if (npoints<1) continue;
+      // cuts on radius      
+      if (npoints==1){
+       if (radius[0]<cradius0||radius[0]>cradius1) continue;
+      }
+      else{
+       if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
+      }
+      //      
+      Double_t delta1=10000,delta2=10000;
+      // cuts on the intersection radius
+      helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
+      if (radius[0]<20&&delta1<1) continue; //intersection at vertex
+      if (radius[0]<10&&delta1<3) continue; //intersection at vertex
+      if (npoints==2){ 
+       helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
+       if (radius[1]<20&&delta2<1) continue;  //intersection at vertex
+       if (radius[1]<10&&delta2<3) continue;  //intersection at vertex 
+      }
+      //
+      Double_t distance1 = TMath::Min(delta1,delta2);
+      if (distance1>cdist1) continue;  // cut on DCA linear approximation
+      //
+      npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
+      helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
+      if (radius[0]<20&&delta1<1) continue; //intersection at vertex
+      if (radius[0]<10&&delta1<3) continue; //intersection at vertex
+      //
+      if (npoints==2){ 
+       helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);   
+       if (radius[1]<20&&delta2<1) continue;  //intersection at vertex
+       if (radius[1]<10&&delta2<3) continue;  //intersection at vertex 
+      }            
+      distance1 = TMath::Min(delta1,delta2);
+      Float_t rkink =0;
+      if (delta1<delta2){
+       rkink = TMath::Sqrt(radius[0]);
+      }
+      else{
+       rkink = TMath::Sqrt(radius[1]);
+      }
+      if (distance1>cdist2) continue;
+      //
+      //
+      AliTPCseed * track1 = (AliTPCseed*)array->At(j);
+      //
+      //
+      Int_t row0 = GetRowNumber(rkink); 
+      if (row0<10)  continue;
+      if (row0>150) continue;
+      //
+      //
+      Float_t dens00=-1,dens01=-1;
+      Float_t dens10=-1,dens11=-1;
+      //
+      Int_t found,foundable,shared;
+      track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
+      if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
+      track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
+      if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
+      //
+      track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
+      if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
+      track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
+      if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
+      //
+      if (dens00<dens10 && dens01<dens11) continue;
+      if (dens00>dens10 && dens01>dens11) continue;
+      if (TMath::Max(dens00,dens10)<0.1)  continue;
+      if (TMath::Max(dens01,dens11)<0.3)  continue;
+      //
+      if (TMath::Min(dens00,dens10)>0.6)  continue;
+      if (TMath::Min(dens01,dens11)>0.6)  continue;
+
+      //
+      AliTPCseed * ktrack0, *ktrack1;
+      if (dens00>dens10){
+       ktrack0 = track0;
+       ktrack1 = track1;
+      }
+      else{
+       ktrack0 = track1;
+       ktrack1 = track0;
+      }
+      if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
+      AliExternalTrackParam paramm(*ktrack0);
+      AliExternalTrackParam paramd(*ktrack1);
+      if (row0>60&&ktrack1->GetReference().X()>90.) new (&paramd) AliExternalTrackParam(ktrack1->GetReference()); 
+      //
+      //
+      kink->SetMother(paramm);
+      kink->SetDaughter(paramd);
+      kink->Update();
+
+      Float_t x[3] = { kink->fXr[0],kink->fXr[1],kink->fXr[2]};
+      Int_t index[4];
+      fParam->Transform0to1(x,index);
+      fParam->Transform1to2(x,index);
+      row0 = GetRowNumber(x[0]); 
+
+      if (kink->fRr<100) continue;
+      if (kink->fRr>240) continue;
+      if (kink->fDist2>cdist3) continue;
+      Float_t dird = kink->fPdr[0]*kink->fXr[0]+kink->fPdr[1]*kink->fXr[1];  // rough direction estimate
+      if (dird<0) continue;
+
+      Float_t dirm = kink->fPm[0]*kink->fXr[0]+kink->fPm[1]*kink->fXr[1];  // rough direction estimate
+      if (dirm<0) continue;
+      Float_t mpt = TMath::Sqrt(kink->fPm[0]*kink->fPm[0]+kink->fPm[1]*kink->fPm[1]);
+      if (mpt<0.2) continue;
+
+
+      Double_t qt   =  TMath::Sin(kink->fAngle[2])*ktrack1->P();
+      if (qt>0.35) continue; 
+      
+      kink->fLab[0] = CookLabel(ktrack0,0.4,0,row0);
+      kink->fLab[1] = CookLabel(ktrack1,0.4,row0,160);
+      if (dens00>dens10){
+       kink->fTPCdensity[0][0] = dens00;
+       kink->fTPCdensity[0][1] = dens01;
+       kink->fTPCdensity[1][0] = dens10;
+       kink->fTPCdensity[1][1] = dens11;
+       kink->fIndex[0] = i;
+       kink->fIndex[1] = j;
+       kink->fZm[0] = zm[i];
+       kink->fZm[1] = zm[j];
+       kink->fFi[0] = fim[i];
+       kink->fFi[1] = fim[j];
+      }
+      else{
+       kink->fTPCdensity[0][0] = dens10;
+       kink->fTPCdensity[0][1] = dens11;
+       kink->fTPCdensity[1][0] = dens00;
+       kink->fTPCdensity[1][1] = dens01;
+       kink->fIndex[0] = j;
+       kink->fIndex[1] = i;
+       kink->fZm[0] = zm[j];
+       kink->fZm[1] = zm[i];
+       kink->fFi[0] = fim[j];
+       kink->fFi[1] = fim[i];
+
+      }
+      if (kink->GetTPCDensityFactor()<0.8) continue;
+      if ((2-kink->GetTPCDensityFactor())*kink->fDist2 >0.25) continue;
+      if (kink->fAngle[2]*ktrack0->P()<0.003) continue; //too small angle
+      if (kink->fAngle[2]>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
+      if (kink->fAngle[2]>0.2&&kink->fTPCdensity[0][1]) continue;
+      if (kink->fAngle[2]<0.02) continue;
+
+
+      //
+      Int_t drow = Int_t(2./TMath::Tan(0.3+kink->fAngle[2]));  // overlap region defined
+      kink->fRow0 = row0;
+      Float_t shapesum =0;
+      Float_t sum = 0;
+      for ( Int_t row = row0-drow; row<row0+drow;row++){
+       if (row<0) continue;
+       if (row>155) continue;
+       if (ktrack0->fClusterPointer[row]){
+         AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
+         shapesum+=point->GetSigmaY()+point->GetSigmaZ();
+         sum++;
+       }
+       if (ktrack1->fClusterPointer[row]){
+         AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
+         shapesum+=point->GetSigmaY()+point->GetSigmaZ();
+         sum++;
+       }       
+      }
+      if (sum<4){
+       kink->fShapeFactor=-1;
+      }
+      else{
+       kink->fShapeFactor = shapesum/sum;
+      }
+      
+      //      esd->AddKink(kink);
+      kinks->AddLast(kink);
+      kink = new AliESDkink;
+      ncandidates++;
+    }
+  }
+  Int_t nkinks = kinks->GetEntriesFast();
+  Float_t *quality     = new Float_t[nkinks];
+  Int_t   * indexes = new Int_t[nkinks];
+  for (Int_t i=0;i<nkinks;i++){
+    quality[i] =100000;
+    AliESDkink *kink = (AliESDkink*)kinks->At(i);
+    if (kink) quality[i] = kink->fDist2*(2.-kink->GetTPCDensityFactor());
+  }
+  TMath::Sort(nkinks,quality,indexes,kFALSE);
+  
+  for (Int_t i=0;i<nkinks;i++){
+    AliESDkink * kink = (AliESDkink*) kinks->At(indexes[i]);
+    Int_t index0 = kink->fIndex[0];
+    Int_t index1 = kink->fIndex[1];
+    kink->fMultiple[0] = usage[index0];
+    kink->fMultiple[1] = usage[index1];
+    if (kink->fMultiple[0]+kink->fMultiple[1]>2) continue;
+    if (kink->fMultiple[0]+kink->fMultiple[1]>0 && quality[indexes[i]]>0.2) continue;
+
+    Int_t index = esd->AddKink(kink);
+    AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
+    AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
+    ktrack0->fKinkIndexes[usage[index0]] = -(index+1);
+    ktrack1->fKinkIndexes[usage[index1]] =  (index+1);
+    usage[index0]++;
+    usage[index1]++;
+  }
+  delete [] quality;
+  delete [] indexes;
+
+  delete[] fim;
+  delete[] zm;
+  delete [] usage;
+  delete[] alpha;
+  delete[] nclusters;
+  delete[] sign;
+  delete[] helixes;
+  kinks->Delete();
+  delete kinks;
+
+  printf("Ncandidates=\t%d\t%d\t%d\n",ncandidates,ntracks,nall);
+  timer.Print();
+}
+
+
+
+void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
+  //
+  // update Kink quality information for mother after back propagation
+  //
+  if (seed->GetKinkIndex(0)>=0) return; 
+  for (Int_t ikink=0;ikink<3;ikink++){
+    Int_t index = seed->GetKinkIndex(ikink);
+    if (index>=0) break;
+    index = TMath::Abs(index)-1;
+    AliESDkink * kink = fEvent->GetKink(index);
+    kink->fTPCdensity2[0][0]=-1;
+    kink->fTPCdensity2[0][1]=-1;
+    //
+    Int_t row0 = kink->fRow0 - Int_t( 2./ (0.05+kink->fAngle[2]));
+    if (row0<15) row0=15;
+    //
+    Int_t row1 = kink->fRow0 + Int_t( 2./ (0.05+kink->fAngle[2]));
+    if (row1>145) row1=145;
+    //
+    Int_t found,foundable,shared;
+    seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
+    if (foundable>5)   kink->fTPCdensity2[0][0] = Float_t(found)/Float_t(foundable);
+    seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
+    if (foundable>5)   kink->fTPCdensity2[0][1] = Float_t(found)/Float_t(foundable);
+    if (kink->fTPCdensity2[0][1]>0.5) kink->fStatus=-100000;
+    if (kink->fTPCdensity2[0][0]<0.5) kink->fStatus=-100000;
+    if (kink->fTPCdensity2[0][0]-kink->fTPCdensity2[0][1]<0.2) kink->fStatus=-100000;
+
+    if (kink->fDist2>1) kink->fStatus=-10000;
+  }
+    
+}
+
 Int_t  AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed, Float_t th)
 {
   //
@@ -4080,8 +4605,12 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() {
       }
     } 
   }
-  RemoveDouble(fSeeds,0.2,0.6,11);
-  RemoveUsed(fSeeds,0.5,0.5,6);
+  //
+  RemoveUsed2(fSeeds,0.85,0.85,0);
+  FindKinks(fSeeds,fEvent);
+  RemoveUsed2(fSeeds,0.5,0.4,30);
+  //RemoveDouble(fSeeds,0.2,0.6,11);
+  //  RemoveUsed(fSeeds,0.5,0.5,6);
 
   //
   Int_t nseed=fSeeds->GetEntriesFast();
@@ -4458,6 +4987,7 @@ void  AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rla
     for (Int_t i=0; i<nseed; i++) {
       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;       
       if (!pt) continue;
+      if (nr==80) pt->UpdateReference();
       if (!pt->IsActive()) continue;
       //      if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
       if (pt->fRelativeSector>17) {
@@ -4542,7 +5072,7 @@ Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
   Int_t nseed= arr->GetEntriesFast();
   for (Int_t i=0;i<nseed;i++){
     AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
-    if (pt) { 
+    if (pt&& pt->GetKinkIndex(0)<=0) { 
       //AliTPCseed *pt2 = new AliTPCseed(*pt);
       fSectors = fInnerSec;
       //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
@@ -4552,7 +5082,14 @@ Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
       //       Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
       //       FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
       //}
-    }      
+    }
+    if (pt&& pt->GetKinkIndex(0)>0) {
+      AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
+      pt->fFirstPoint = kink->fRow0;
+      fSectors = fInnerSec;
+      FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);  
+    }
+    
   }
   return 0;
 }
@@ -4692,8 +5229,6 @@ Float_t  AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
 }
 
 
-
-
 //__________________________________________________________________________
 void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const {
   //--------------------------------------------------------------------
@@ -4777,6 +5312,90 @@ void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const {
 }
 
 
+//__________________________________________________________________________
+Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
+  //--------------------------------------------------------------------
+  //This function "cooks" a track label. If label<0, this track is fake.
+  //--------------------------------------------------------------------
+  Int_t noc=t->GetNumberOfClusters();
+  if (noc<10){
+    //printf("\nnot founded prolongation\n\n\n");
+    //t->Dump();
+    return -1;
+  }
+  Int_t lb[160];
+  Int_t mx[160];
+  AliTPCclusterMI *clusters[160];
+  //
+  for (Int_t i=0;i<160;i++) {
+    clusters[i]=0;
+    lb[i]=mx[i]=0;
+  }
+
+  Int_t i;
+  Int_t current=0;
+  for (i=0; i<160 && current<noc; i++) {
+    if (i<first) continue;
+    if (i>last)  continue;
+     Int_t index=t->GetClusterIndex2(i);
+     if (index<=0) continue; 
+     if (index&0x8000) continue;
+     //     
+     //clusters[current]=GetClusterMI(index);
+     if (t->fClusterPointer[i]){
+       clusters[current]=t->fClusterPointer[i];     
+       current++;
+     }
+  }
+  noc = current;
+  if (noc<5) return -1;
+  Int_t lab=123456789;
+  for (i=0; i<noc; i++) {
+    AliTPCclusterMI *c=clusters[i];
+    if (!c) continue;
+    lab=TMath::Abs(c->GetLabel(0));
+    Int_t j;
+    for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
+    lb[j]=lab;
+    (mx[j])++;
+  }
+
+  Int_t max=0;
+  for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
+    
+  for (i=0; i<noc; i++) {
+    AliTPCclusterMI *c=clusters[i]; 
+    if (!c) continue;
+    if (TMath::Abs(c->GetLabel(1)) == lab ||
+        TMath::Abs(c->GetLabel(2)) == lab ) max++;
+  }
+
+  if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
+
+  else {
+     Int_t tail=Int_t(0.10*noc);
+     max=0;
+     Int_t ind=0;
+     for (i=1; i<=160&&ind<tail; i++) {
+       //       AliTPCclusterMI *c=clusters[noc-i];
+       AliTPCclusterMI *c=clusters[i];
+       if (!c) continue;
+       if (lab == TMath::Abs(c->GetLabel(0)) ||
+           lab == TMath::Abs(c->GetLabel(1)) ||
+           lab == TMath::Abs(c->GetLabel(2))) max++;
+       ind++;
+     }
+     if (max < Int_t(0.5*tail)) lab=-lab;
+  }
+
+  //  t->SetLabel(lab);
+  return lab;
+  //  delete[] lb;
+  //delete[] mx;
+  //delete[] clusters;
+}
+
+
 Int_t  AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const 
 {
   //return pad row number for this x
@@ -4984,7 +5603,7 @@ AliTPCseed::AliTPCseed():AliTPCtrack(){
   fRemoval =0; 
   for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
-
+  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
   fPoints = 0;
   fEPoints = 0;
   fNFoundable =0;
@@ -5014,6 +5633,8 @@ AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
   //  fTrackPoints =0;
   fRemoval =0;
   fSort =0;
+  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=t.GetKinkIndex(i);
+
   for (Int_t i=0;i<160;i++) {
     fClusterPointer[i] = 0;
     Int_t index = t.GetClusterIndex(i);
@@ -5043,7 +5664,7 @@ AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
     Int_t index = t.GetClusterIndex(i);
     SetClusterIndex2(i,index);
   }
-  
+  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
   fPoints = 0;
   fEPoints = 0;
   fNFoundable =0; 
@@ -5062,15 +5683,17 @@ AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
 
 }
 
-AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15], 
-                                       Double_t xr, Double_t alpha):      
+
+AliTPCseed::AliTPCseed(UInt_t index,  const Double_t xx[5], const Double_t cc[15], 
+                                        Double_t xr, Double_t alpha):      
   AliTPCtrack(index, xx, cc, xr, alpha) {
-  //
+   //
   //
   //constructor
   fRow =0;
   for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
+  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
   fPoints = 0;
   fEPoints = 0;
   fNFoundable =0;
index 4a13d1e75252d6aefa977726d05d3a5e002b7095..be516cebb29e6a04ad548266197aa9cc476233b7 100644 (file)
@@ -35,6 +35,8 @@ class AliTPCseed : public AliTPCtrack {
      AliTPCseed(const AliTPCtrack &t);
      AliTPCseed(const AliTPCseed &s);
      AliTPCseed(const AliKalmanTrack &t, Double_t a);
+     AliTPCseed(UInt_t index, const Double_t xx[5], 
+                const Double_t cc[15], Double_t xr, Double_t alpha);     
      Int_t Compare(const TObject *o) const;
      void Reset(Bool_t all = kTRUE);
      Int_t GetProlongation(Double_t xr, Double_t &y, Double_t & z) const;
@@ -42,9 +44,6 @@ class AliTPCseed : public AliTPCtrack {
      virtual Int_t Update(const AliTPCclusterMI* c, Double_t chi2, UInt_t i);
      AliTPCTrackerPoint * GetTrackPoint(Int_t i);
      void RebuildSeed(); // rebuild seed to be ready for storing
-     AliTPCseed(UInt_t index, const Double_t xx[5], 
-                const Double_t cc[15], Double_t xr, Double_t alpha);
-
      Double_t GetDensityFirst(Int_t n);
      Double_t GetSigma2C() const {return fC44;};
      void GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2);
@@ -136,6 +135,9 @@ public:
   void WriteTracks(TTree * tree);  
   void DeleteSeeds();
   void SetDebug(Int_t debug){ fDebug = debug;}
+  void FindKinks(TObjArray * array, AliESD * esd);
+  void UpdateKinkQualityM(AliTPCseed * seed);
+
    Int_t ReadSeeds(const TFile *in);
    TObjArray * GetSeeds(){return fSeeds;}
    //   
@@ -143,6 +145,7 @@ public:
    AliTPCclusterMI *GetClusterMI(Int_t index) const;
    Int_t Clusters2Tracks();
    virtual void  CookLabel(AliTPCseed *t,Float_t wrong) const; 
+   virtual Int_t   CookLabel(AliTPCseed *t,Float_t wrong, Int_t first,Int_t last ) const; 
    
    void RotateToLocal(AliTPCseed *seed);
   
@@ -260,6 +263,7 @@ private:
    void  SignShared(TObjArray * arr);
 
    void  RemoveUsed(TObjArray * arr, Float_t factor1, Float_t factor2,  Int_t removalindex);
+   void  RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal);
    void  RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2,  Int_t removalindex);
 
    void  StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const;
index e59a5230e4388778a6c1805c3a18b3a45e5c62a9..167a34f87674fd5c5a9998e5c0edc6c887564599 100644 (file)
@@ -32,7 +32,7 @@
 #pragma link C++ class AliTPCtrack+;
 #pragma link C++ class AliTPCtracker+;
 
-#pragma link C++ class AliHelix+;
+//#pragma link C++ class AliHelix+;
 #pragma link C++ class AliTPCpolyTrack+;
 #pragma link C++ class AliTPCseed+;           // defined in AliTPCtrackerMI.h
 #pragma link C++ class AliTPCtrackerMI+;
index 71654b469cbc8c14536b75653e6ef0230d081002..a61fe3a360634b9fc1ef00108bc7731658471bc3 100644 (file)
@@ -4,7 +4,7 @@ SRCS:=  AliTPCcluster.cxx \
         AliClustersArray.cxx AliTPCClustersArray.cxx \
        AliTPCclusterer.cxx AliTPCclustererMI.cxx \
        AliTPCtrack.cxx AliTPCtracker.cxx \
-       AliHelix.cxx AliTPCpolyTrack.cxx AliTPCtrackerMI.cxx \
+       AliTPCpolyTrack.cxx AliTPCtrackerMI.cxx \
        AliTPCPid.cxx AliTPCtrackPid.cxx AliTPCpidESD.cxx \
        AliTPCReconstructor.cxx