]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliESDtrack.cxx
Initialization of fAmpThreshold is added to ctors.
[u/mrichter/AliRoot.git] / STEER / AliESDtrack.cxx
index 4ed8ca6ecbdbd6407b3633cb5027c6ea51078ea0..44669f7b2c291cb0548965c035f244f5ac39db64 100644 (file)
 //   This is the class to deal with during the phisics analysis of data
 //      Origin: Iouri Belikov, CERN
 //      e-mail: Jouri.Belikov@cern.ch
+//
+//
+//
+//  What do you need to know before starting analysis
+//  (by Marian Ivanov: marian.ivanov@cern.ch)
+//
+//
+//   AliESDtrack:
+//   1.  What is the AliESDtrack
+//   2.  What informations do we store
+//   3.  How to use the information for analysis
+//   
+//
+//   1.AliESDtrack is the container of the information about the track/particle
+//     reconstructed during Barrel Tracking.
+//     The track information is propagated from one tracking detector to 
+//     other using the functionality of AliESDtrack - Current parameters.  
+//
+//     No global fit model is used.
+//     Barrel tracking use Kalman filtering technique, it gives optimal local 
+//     track parameters at given point under certian assumptions.
+//             
+//     Kalman filter take into account additional effect which are 
+//     difficult to handle using global fit.
+//     Effects:
+//        a.) Multiple scattering
+//        b.) Energy loss
+//        c.) Non homogenous magnetic field
+//
+//     In general case, following barrel detectors are contributing to 
+//     the Kalman track information:
+//         a. TPC
+//         b. ITS
+//         c. TRD
+//
+//      In general 3 reconstruction itteration are performed:
+//         1. Find tracks   - sequence TPC->ITS
+//         2. PropagateBack - sequence ITS->TPC->TRD -> Outer PID detectors
+//         3. Refit invward - sequence TRD->TPC->ITS
+//      The current tracks are updated after each detector (see bellow).
+//      In specical cases a track  sanpshots are stored.   
+// 
+//
+//      For some type of analysis (+visualization) track local parameters at 
+//      different position are neccesary. A snapshots during the track 
+//      propagation are created.
+//      (See AliExternalTrackParam class for desctiption of variables and 
+//      functionality)
+//      Snapshots:
+//      a. Current parameters - class itself (AliExternalTrackParam)
+//         Contributors: general case TRD->TPC->ITS
+//         Preferable usage:  Decission  - primary or secondary track
+//         NOTICE - By default the track parameters are stored at the DCA point
+//                  to the primary vertex. optimal for primary tracks, 
+//                  far from optimal for secondary tracks.
+//      b. Constrained parameters - Kalman information updated with 
+//         the Primary vertex information 
+//         Contributors: general case TRD->TPC->ITS
+//         Preferable usage: Use only for tracks selected as primary
+//         NOTICE - not real constrain - taken as additional measurement 
+//         with corresponding error
+//         Function:  
+//       const AliExternalTrackParam *GetConstrainedParam() const {return fCp;}
+//      c. Inner parameters -  Track parameters at inner wall of the TPC 
+//         Contributors: general case TRD->TPC
+//         function:
+//           const AliExternalTrackParam *GetInnerParam() const { return fIp;}
+//
+//      d. TPCinnerparam  - contributors - TPC only
+//         Contributors:  TPC
+//         Preferable usage: Requested for HBT study 
+//                         (smaller correlations as using also ITS information)
+//         NOTICE - the track parameters are propagated to the DCA to  
+//         to primary vertex
+//         Optimal for primary, far from optimal for secondary tracks
+//         Function:
+//    const AliExternalTrackParam *GetTPCInnerParam() const {return fTPCInner;}
+//     
+//      e. Outer parameters - 
+//           Contributors-  general case - ITS-> TPC -> TRD
+//           The last point - Outer parameters radius is determined
+//           e.a) Local inclination angle bigger than threshold - 
+//                Low momenta tracks 
+//           e.a) Catastrofic energy losss in material
+//           e.b) Not further improvement (no space points)
+//           Usage:             
+//              a.) Tracking: Starting parameter for Refit inward 
+//              b.) Visualization
+//              c.) QA
+//         NOTICE: Should be not used for the physic analysis
+//         Function:
+//            const AliExternalTrackParam *GetOuterParam() const { return fOp;}
+//
 //-----------------------------------------------------------------
 
 #include <TMath.h>
 #include "AliESDVertex.h"
 #include "AliESDtrack.h"
 #include "AliKalmanTrack.h"
+#include "AliVTrack.h"
 #include "AliLog.h"
 #include "AliTrackPointArray.h"
+#include "TPolyMarker3D.h"
 
 ClassImp(AliESDtrack)
 
@@ -75,14 +170,17 @@ AliESDtrack::AliESDtrack() :
   fTPCLabel(0),
   fTRDLabel(0),
   fTOFCalChannel(0),
-  fTOFindex(0),
+  fTOFindex(-1),
   fHMPIDqn(0),
-  fHMPIDcluIdx(0),
+  fHMPIDcluIdx(-1),
   fEMCALindex(kEMCALNoMatch),
   fHMPIDtrkTheta(0),
   fHMPIDtrkPhi(0),
   fHMPIDsignal(0),
   fTrackLength(0),
+  fdTPC(0),fzTPC(0),
+  fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
+  fCchi2TPC(0),
   fD(0),fZ(0),
   fCdd(0),fCdz(0),fCzz(0),
   fCchi2(0),
@@ -91,6 +189,7 @@ AliESDtrack::AliESDtrack() :
   fTRDchi2(0),
   fTOFchi2(0),
   fHMPIDchi2(0),
+  fGlobalChi2(0),
   fITSsignal(0),
   fTPCsignal(0),
   fTPCsignalS(0),
@@ -112,12 +211,15 @@ AliESDtrack::AliESDtrack() :
   fITSClusterMap(0),
   fTRDncls(0),
   fTRDncls0(0),
-  fTRDpidQuality(0)
+  fTRDpidQuality(0),
+  fTRDnSlices(0),
+  fTRDslices(0x0)
+  
 {
   //
   // The default ESD constructor 
   //
-  Int_t i, j;
+  Int_t i;
   for (i=0; i<AliPID::kSPECIES; i++) {
     fTrackTime[i]=0.;
     fR[i]=0.;
@@ -130,15 +232,13 @@ AliESDtrack::AliESDtrack() :
   
   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
   for (i=0; i<3; i++)   { fV0Indexes[i]=0;}
-  for (i=0;i<kNPlane;i++) {
-    for (j=0;j<kNSlice;j++) {
-      fTRDsignals[i][j]=0.; 
-    }
+  for (i=0;i<kTRDnPlanes;i++) {
     fTRDTimBin[i]=0;
   }
   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
+  for (i=0;i<12;i++) {fITSModule[i]=-1;}
 }
 
 //_______________________________________________________________________
@@ -166,6 +266,9 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
   fHMPIDtrkPhi(track.fHMPIDtrkPhi),
   fHMPIDsignal(track.fHMPIDsignal),
   fTrackLength(track.fTrackLength),
+  fdTPC(track.fdTPC),fzTPC(track.fzTPC),
+  fCddTPC(track.fCddTPC),fCdzTPC(track.fCdzTPC),fCzzTPC(track.fCzzTPC),
+  fCchi2TPC(track.fCchi2TPC),
   fD(track.fD),fZ(track.fZ),
   fCdd(track.fCdd),fCdz(track.fCdz),fCzz(track.fCzz),
   fCchi2(track.fCchi2),
@@ -174,6 +277,7 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
   fTRDchi2(track.fTRDchi2),
   fTOFchi2(track.fTOFchi2),
   fHMPIDchi2(track.fHMPIDchi2),
+  fGlobalChi2(track.fGlobalChi2),
   fITSsignal(track.fITSsignal),
   fTPCsignal(track.fTPCsignal),
   fTPCsignalS(track.fTPCsignalS),
@@ -195,7 +299,9 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
   fITSClusterMap(track.fITSClusterMap),
   fTRDncls(track.fTRDncls),
   fTRDncls0(track.fTRDncls0),
-  fTRDpidQuality(track.fTRDpidQuality)
+  fTRDpidQuality(track.fTRDpidQuality),
+  fTRDnSlices(track.fTRDnSlices),
+  fTRDslices(0x0)
 {
   //
   //copy constructor
@@ -210,16 +316,20 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
   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];}
   //
-  for (Int_t i=0;i<kNPlane;i++) {
-    for (Int_t j=0;j<kNSlice;j++) {
-      fTRDsignals[i][j]=track.fTRDsignals[i][j]; 
-    }
+  for (Int_t i=0;i<kTRDnPlanes;i++) {
     fTRDTimBin[i]=track.fTRDTimBin[i];
   }
+
+  if (fTRDnSlices) {
+    fTRDslices=new Double32_t[fTRDnSlices];
+    for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=track.fTRDslices[i];
+  }
+
   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i]=track.fTRDr[i]; 
   for (Int_t i=0;i<AliPID::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<12;i++) fITSModule[i]=track.fITSModule[i];
   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i]=track.fHMPIDr[i];
 
   if (track.fCp) fCp=new AliExternalTrackParam(*track.fCp);
@@ -230,6 +340,124 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
   if (track.fFriendTrack) fFriendTrack=new AliESDfriendTrack(*(track.fFriendTrack));
 }
 
+//_______________________________________________________________________
+AliESDtrack::AliESDtrack(const AliVTrack *track) : 
+  AliExternalTrackParam(track),
+  fCp(0),
+  fIp(0),
+  fTPCInner(0),
+  fOp(0),
+  fFriendTrack(0),
+  fTPCClusterMap(159),//number of padrows
+  fTPCSharedMap(159),//number of padrows
+  fFlags(0),
+  fID(),
+  fLabel(0),
+  fITSLabel(0),
+  fTPCLabel(0),
+  fTRDLabel(0),
+  fTOFCalChannel(0),
+  fTOFindex(-1),
+  fHMPIDqn(0),
+  fHMPIDcluIdx(-1),
+  fEMCALindex(kEMCALNoMatch),
+  fHMPIDtrkTheta(0),
+  fHMPIDtrkPhi(0),
+  fHMPIDsignal(0),
+  fTrackLength(0),
+  fdTPC(0),fzTPC(0),
+  fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
+  fCchi2TPC(0),
+  fD(0),fZ(0),
+  fCdd(0),fCdz(0),fCzz(0),
+  fCchi2(0),
+  fITSchi2(0),
+  fTPCchi2(0),
+  fTRDchi2(0),
+  fTOFchi2(0),
+  fHMPIDchi2(0),
+  fGlobalChi2(0),
+  fITSsignal(0),
+  fTPCsignal(0),
+  fTPCsignalS(0),
+  fTRDsignal(0),
+  fTRDQuality(0),
+  fTRDBudget(0),
+  fTOFsignal(0),
+  fTOFsignalToT(0),
+  fTOFsignalRaw(0),
+  fTOFsignalDz(0),
+  fHMPIDtrkX(0),
+  fHMPIDtrkY(0),
+  fHMPIDmipX(0),
+  fHMPIDmipY(0),
+  fTPCncls(0),
+  fTPCnclsF(0),
+  fTPCsignalN(0),
+  fITSncls(0),
+  fITSClusterMap(0),
+  fTRDncls(0),
+  fTRDncls0(0),
+  fTRDpidQuality(0),
+  fTRDnSlices(0),
+  fTRDslices(0x0)
+{
+  //
+  // ESD track from AliVTrack.
+  // This is not a copy constructor !
+  //
+
+  if (track->InheritsFrom("AliExternalTrackParam")) {
+     AliError("This is not a copy constructor. Use AliESDtrack(const AliESDtrack &) !");
+     AliWarning("Calling the default constructor...");
+     AliESDtrack();
+     return;
+  }
+
+  // Reset all the arrays
+  Int_t i;
+  for (i=0; i<AliPID::kSPECIES; i++) {
+    fTrackTime[i]=0.;
+    fR[i]=0.;
+    fITSr[i]=0.;
+    fTPCr[i]=0.;
+    fTRDr[i]=0.;
+    fTOFr[i]=0.;
+    fHMPIDr[i]=0.;
+  }
+  
+  for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
+  for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
+  for (i=0;i<kTRDnPlanes;i++) {
+    fTRDTimBin[i]=0;
+  }
+  for (i=0;i<4;i++) {fTPCPoints[i]=0;}
+  for (i=0;i<3;i++) {fTOFLabel[i]=0;}
+  for (i=0;i<10;i++) {fTOFInfo[i]=0;}
+  for (i=0;i<12;i++) {fITSModule[i]=-1;}
+
+  // Set the ID
+  SetID(track->GetID());
+
+  // Set ITS cluster map
+  fITSClusterMap=track->GetITSClusterMap();
+
+  fITSncls=0;
+  for(i=0; i<6; i++) {
+    if(HasPointOnITSLayer(i)) fITSncls++;
+  }
+
+  // Set the combined PID
+  const Double_t *pid = track->PID();
+  if(pid){
+    for (i=0; i<AliPID::kSPECIES; i++) fR[i]=pid[i];
+  }
+  // AliESD track label
+  SetLabel(track->GetLabel());
+  // Set the status
+  SetStatus(track->GetStatus());
+}
+
 //_______________________________________________________________________
 AliESDtrack::AliESDtrack(TParticle * part) : 
   AliExternalTrackParam(),
@@ -247,14 +475,17 @@ AliESDtrack::AliESDtrack(TParticle * part) :
   fTPCLabel(0),
   fTRDLabel(0),
   fTOFCalChannel(0),
-  fTOFindex(0),
+  fTOFindex(-1),
   fHMPIDqn(0),
-  fHMPIDcluIdx(0),
+  fHMPIDcluIdx(-1),
   fEMCALindex(kEMCALNoMatch),
   fHMPIDtrkTheta(0),
   fHMPIDtrkPhi(0),
   fHMPIDsignal(0),
   fTrackLength(0),
+  fdTPC(0),fzTPC(0),
+  fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
+  fCchi2TPC(0),
   fD(0),fZ(0),
   fCdd(0),fCdz(0),fCzz(0),
   fCchi2(0),
@@ -263,6 +494,7 @@ AliESDtrack::AliESDtrack(TParticle * part) :
   fTRDchi2(0),
   fTOFchi2(0),
   fHMPIDchi2(0),
+  fGlobalChi2(0),
   fITSsignal(0),
   fTPCsignal(0),
   fTPCsignalS(0),
@@ -284,14 +516,16 @@ AliESDtrack::AliESDtrack(TParticle * part) :
   fITSClusterMap(0),
   fTRDncls(0),
   fTRDncls0(0),
-  fTRDpidQuality(0)
+  fTRDpidQuality(0),
+  fTRDnSlices(0),
+  fTRDslices(0x0)
 {
   //
   // ESD track from TParticle
   //
 
   // Reset all the arrays
-  Int_t i, j;
+  Int_t i;
   for (i=0; i<AliPID::kSPECIES; i++) {
     fTrackTime[i]=0.;
     fR[i]=0.;
@@ -304,15 +538,13 @@ AliESDtrack::AliESDtrack(TParticle * part) :
   
   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
-  for (i=0;i<kNPlane;i++) {
-    for (j=0;j<kNSlice;j++) {
-      fTRDsignals[i][j]=0.; 
-    }
+  for (i=0;i<kTRDnPlanes;i++) {
     fTRDTimBin[i]=0;
   }
   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
+  for (i=0;i<12;i++) {fITSModule[i]=-1;}
 
   // Calculate the AliExternalTrackParam content
 
@@ -414,8 +646,204 @@ AliESDtrack::~AliESDtrack(){
   delete fOp;
   delete fCp; 
   delete fFriendTrack;
+  if(fTRDnSlices)
+    delete[] fTRDslices;
 }
 
+AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
+  
+
+  if(&source == this) return *this;
+  AliExternalTrackParam::operator=(source);
+
+  
+  if(source.fCp){
+    // we have the trackparam: assign or copy construct
+    if(fCp)*fCp = *source.fCp;
+    else fCp = new AliExternalTrackParam(*source.fCp);
+  }
+  else{
+    // no track param delete the old one
+    if(fCp)delete fCp;
+    fCp = 0;
+  }
+
+  if(source.fIp){
+    // we have the trackparam: assign or copy construct
+    if(fIp)*fIp = *source.fIp;
+    else fIp = new AliExternalTrackParam(*source.fIp);
+  }
+  else{
+    // no track param delete the old one
+    if(fIp)delete fIp;
+    fIp = 0;
+  }
+
+
+  if(source.fTPCInner){
+    // we have the trackparam: assign or copy construct
+    if(fTPCInner) *fTPCInner = *source.fTPCInner;
+    else fTPCInner = new AliExternalTrackParam(*source.fTPCInner);
+  }
+  else{
+    // no track param delete the old one
+    if(fTPCInner)delete fTPCInner;
+    fTPCInner = 0;
+  }
+
+
+  if(source.fOp){
+    // we have the trackparam: assign or copy construct
+    if(fOp) *fOp = *source.fOp;
+    else fOp = new AliExternalTrackParam(*source.fOp);
+  }
+  else{
+    // no track param delete the old one
+    if(fOp)delete fOp;
+    fOp = 0;
+  }
+
+  // copy also the friend track 
+  // use copy constructor
+  if(source.fFriendTrack){
+    // we have the trackparam: assign or copy construct
+    delete fFriendTrack; fFriendTrack=new AliESDfriendTrack(*source.fFriendTrack);
+  }
+  else{
+    // no track param delete the old one
+    delete fFriendTrack; fFriendTrack= 0;
+  }
+
+  fTPCClusterMap = source.fTPCClusterMap; 
+  fTPCSharedMap  = source.fTPCSharedMap;  
+  // the simple stuff
+  fFlags    = source.fFlags; 
+  fID       = source.fID;             
+  fLabel    = source.fLabel;
+  fITSLabel = source.fITSLabel;
+  for(int i = 0; i< 12;++i){
+    fITSModule[i] = source.fITSModule[i];
+  }
+  fTPCLabel = source.fTPCLabel; 
+  fTRDLabel = source.fTRDLabel;
+  for(int i = 0; i< 3;++i){
+    fTOFLabel[i] = source.fTOFLabel[i];    
+  }
+  fTOFCalChannel = source.fTOFCalChannel;
+  fTOFindex      = source.fTOFindex;
+  fHMPIDqn       = source.fHMPIDqn;
+  fHMPIDcluIdx   = source.fHMPIDcluIdx; 
+  fEMCALindex    = source.fEMCALindex;
+
+  for(int i = 0; i< 3;++i){
+    fKinkIndexes[i] = source.fKinkIndexes[i]; 
+    fV0Indexes[i]   = source.fV0Indexes[i]; 
+  }
+
+  for(int i = 0; i< AliPID::kSPECIES;++i){
+    fR[i]     = source.fR[i];
+    fITSr[i]  = source.fITSr[i];
+    fTPCr[i]  = source.fTPCr[i];
+    fTRDr[i]  = source.fTRDr[i];
+    fTOFr[i]  = source.fTOFr[i];
+    fHMPIDr[i] = source.fHMPIDr[i];
+    fTrackTime[i] = source.fTrackTime[i];  
+  }
+
+  fHMPIDtrkTheta = source.fHMPIDtrkTheta;
+  fHMPIDtrkPhi   = source.fHMPIDtrkPhi;
+  fHMPIDsignal   = source.fHMPIDsignal; 
+
+  
+  fTrackLength   = source. fTrackLength;
+  fdTPC  = source.fdTPC; 
+  fzTPC  = source.fzTPC; 
+  fCddTPC = source.fCddTPC;
+  fCdzTPC = source.fCdzTPC;
+  fCzzTPC = source.fCzzTPC;
+  fCchi2TPC = source.fCchi2TPC;
+
+  fD  = source.fD; 
+  fZ  = source.fZ; 
+  fCdd = source.fCdd;
+  fCdz = source.fCdz;
+  fCzz = source.fCzz;
+  fCchi2     = source.fCchi2;
+
+  fITSchi2   = source.fITSchi2;             
+  fTPCchi2   = source.fTPCchi2;            
+  fTRDchi2   = source.fTRDchi2;      
+  fTOFchi2   = source.fTOFchi2;      
+  fHMPIDchi2 = source.fHMPIDchi2;      
+
+  fGlobalChi2 = source.fGlobalChi2;      
+
+  fITSsignal  = source.fITSsignal;     
+  fTPCsignal  = source.fTPCsignal;     
+  fTPCsignalS = source.fTPCsignalS;    
+  for(int i = 0; i< 4;++i){
+    fTPCPoints[i] = source.fTPCPoints[i];  
+  }
+  fTRDsignal = source.fTRDsignal;
+
+  for(int i = 0;i < kTRDnPlanes;++i){
+    fTRDTimBin[i] = source.fTRDTimBin[i];   
+  }
+
+  if(fTRDnSlices)
+    delete[] fTRDslices;
+  fTRDslices=0;
+  fTRDnSlices=source.fTRDnSlices;
+  if (fTRDnSlices) {
+    fTRDslices=new Double32_t[fTRDnSlices];
+    for(int j = 0;j < fTRDnSlices;++j) fTRDslices[j] = source.fTRDslices[j];
+  }
+
+  fTRDQuality =   source.fTRDQuality;     
+  fTRDBudget  =   source.fTRDBudget;      
+  fTOFsignal  =   source.fTOFsignal;     
+  fTOFsignalToT = source.fTOFsignalToT;   
+  fTOFsignalRaw = source.fTOFsignalRaw;  
+  fTOFsignalDz  = source.fTOFsignalDz;      
+  
+  for(int i = 0;i<10;++i){
+    fTOFInfo[i] = source.fTOFInfo[i];    
+  }
+
+  fHMPIDtrkX = source.fHMPIDtrkX; 
+  fHMPIDtrkY = source.fHMPIDtrkY; 
+  fHMPIDmipX = source.fHMPIDmipX;
+  fHMPIDmipY = source.fHMPIDmipY; 
+
+  fTPCncls    = source.fTPCncls;      
+  fTPCnclsF   = source.fTPCnclsF;     
+  fTPCsignalN = source.fTPCsignalN;   
+
+  fITSncls = source.fITSncls;       
+  fITSClusterMap = source.fITSClusterMap; 
+  fTRDncls   = source.fTRDncls;       
+  fTRDncls0  = source.fTRDncls0;      
+  fTRDpidQuality  = source.fTRDpidQuality; 
+  return *this;
+}
+
+
+
+void AliESDtrack::Copy(TObject &obj) const {
+  
+  // this overwrites the virtual TOBject::Copy()
+  // to allow run time copying without casting
+  // in AliESDEvent
+
+  if(this==&obj)return;
+  AliESDtrack *robj = dynamic_cast<AliESDtrack*>(&obj);
+  if(!robj)return; // not an AliESDtrack
+  *robj = *this;
+
+}
+
+
+
 void AliESDtrack::AddCalibObject(TObject * object){
   //
   // add calib object to the list
@@ -433,14 +861,88 @@ TObject *  AliESDtrack::GetCalibObject(Int_t index){
 }
 
 
+Bool_t AliESDtrack::FillTPCOnlyTrack(AliESDtrack &track){
+  
+  // Fills the information of the TPC-only first reconstruction pass
+  // into the passed ESDtrack object. For consistency fTPCInner is also filled
+  // again
+
+
+
+  // For data produced before r26675
+  // RelateToVertexTPC was not properly called during reco
+  // so you'll have to call it again, before FillTPCOnlyTrack
+  //  Float_t p[2],cov[3];
+  // track->GetImpactParametersTPC(p,cov); 
+  // if(p[0]==0&&p[1]==0) // <- Default values
+  //  track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
+  
+
+  if(!fTPCInner)return kFALSE;
+
+  // fill the TPC track params to the global track parameters
+  track.Set(fTPCInner->GetX(),fTPCInner->GetAlpha(),fTPCInner->GetParameter(),fTPCInner->GetCovariance());
+  track.fD = fdTPC;
+  track.fZ = fzTPC;
+  track.fCdd = fCddTPC;
+  track.fCdz = fCdzTPC;
+  track.fCzz = fCzzTPC;
+
+  // copy the TPCinner parameters
+  if(track.fTPCInner) *track.fTPCInner = *fTPCInner;
+  else track.fTPCInner = new AliExternalTrackParam(*fTPCInner);
+  track.fdTPC   = fdTPC;
+  track.fzTPC   = fzTPC;
+  track.fCddTPC = fCddTPC;
+  track.fCdzTPC = fCdzTPC;
+  track.fCzzTPC = fCzzTPC;
+  track.fCchi2TPC = fCchi2TPC;
+
+
+  // copy all other TPC specific parameters
+
+  // replace label by TPC label
+  track.fLabel    = fTPCLabel;
+  track.fTPCLabel = fTPCLabel;
+
+  track.fTPCchi2 = fTPCchi2; 
+  track.fTPCsignal = fTPCsignal;
+  track.fTPCsignalS = fTPCsignalS;
+  for(int i = 0;i<4;++i)track.fTPCPoints[i] = fTPCPoints[i];
+
+  track.fTPCncls    = fTPCncls;     
+  track.fTPCnclsF   = fTPCnclsF;     
+  track.fTPCsignalN =  fTPCsignalN;
+
+  // PID 
+  for(int i=0;i<AliPID::kSPECIES;++i){
+    track.fTPCr[i] = fTPCr[i];
+    // combined PID is TPC only!
+    track.fR[i] = fTPCr[i];
+  }
+  track.fTPCClusterMap = fTPCClusterMap;
+  track.fTPCSharedMap = fTPCSharedMap;
+
+
+  // reset the flags
+  track.fFlags = kTPCin;
+  track.fID    = fID;
+
+  track.fFlags |= fFlags & kTPCpid; //copy the TPCpid status flag
+  for (Int_t i=0;i<3;i++) track.fKinkIndexes[i] = fKinkIndexes[i];
+  
+  return kTRUE;
+    
+}
+
 //_______________________________________________________________________
 void AliESDtrack::MakeMiniESDtrack(){
   // Resets everything except
   // fFlags: Reconstruction status flags 
   // fLabel: Track label
   // fID:  Unique ID of the track
-  // fD: Impact parameter in XY-plane
-  // fZ: Impact parameter in Z 
+  // Impact parameter information
   // fR[AliPID::kSPECIES]: combined "detector response probability"
   // Running track parameters in the base class (AliExternalTrackParam)
   
@@ -450,7 +952,6 @@ void AliESDtrack::MakeMiniESDtrack(){
 
   // Reset track parameters constrained to the primary vertex
   delete fCp;fCp = 0;
-  fCchi2 = 0;
 
   // Reset track parameters at the inner wall of TPC
   delete fIp;fIp = 0;
@@ -487,21 +988,22 @@ void AliESDtrack::MakeMiniESDtrack(){
   fTRDncls = 0;       
   fTRDncls0 = 0;       
   fTRDsignal = 0;      
-  for (Int_t i=0;i<kNPlane;i++) {
-    for (Int_t j=0;j<kNSlice;j++) {
-      fTRDsignals[i][j] = 0; 
-    }
+  for (Int_t i=0;i<kTRDnPlanes;i++) {
     fTRDTimBin[i]  = 0;
   }
   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0; 
   fTRDLabel = 0;       
   fTRDQuality  = 0;
   fTRDpidQuality = 0;
+  if(fTRDnSlices)
+    delete[] fTRDslices;
+  fTRDslices=0x0;
+  fTRDnSlices=0;
   fTRDBudget  = 0;
 
   // Reset TOF related track information
   fTOFchi2 = 0;        
-  fTOFindex = 0;       
+  fTOFindex = -1;       
   fTOFsignal = 0;      
   fTOFCalChannel = 0;
   fTOFsignalToT = 0;
@@ -514,7 +1016,7 @@ void AliESDtrack::MakeMiniESDtrack(){
   // Reset HMPID related track information
   fHMPIDchi2 = 0;     
   fHMPIDqn = 0;     
-  fHMPIDcluIdx = 0;     
+  fHMPIDcluIdx = -1;     
   fHMPIDsignal = 0;     
   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
   fHMPIDtrkTheta = 0;     
@@ -525,6 +1027,9 @@ void AliESDtrack::MakeMiniESDtrack(){
   fHMPIDmipY = 0;
   fEMCALindex = kEMCALNoMatch;
 
+  // reset global track chi2
+  fGlobalChi2 = 0;
+
   delete fFriendTrack; fFriendTrack = 0;
 } 
 //_______________________________________________________________________
@@ -599,6 +1104,7 @@ Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
   switch (flags) {
     
   case kITSin: case kITSout: case kITSrefit:
+    fITSClusterMap=0;
     fITSncls=t->GetNumberOfClusters();
     index=fFriendTrack->GetITSindices(); 
     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
@@ -611,6 +1117,12 @@ Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
     fITSchi2=t->GetChi2();
     fITSsignal=t->GetPIDsignal();
     fITSLabel = t->GetLabel();
+    // keep in fOp the parameters outside ITS for ITS stand-alone tracks 
+    if (flags==kITSout) { 
+      if (!fOp) fOp=new AliExternalTrackParam(*t);
+      else 
+        fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
+    }      
     break;
     
   case kTPCin: case kTPCrefit:
@@ -683,12 +1195,15 @@ Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
     fTPCsignal=t->GetPIDsignal();
     break;
 
-  case kTRDout: case kTRDin: case kTRDrefit:
-    index=fFriendTrack->GetTRDindices();
+  case kTRDin: case kTRDrefit:
+    break;
+  case kTRDout:
+    index     = fFriendTrack->GetTRDindices();
     fTRDLabel = t->GetLabel(); 
-    fTRDncls=t->GetNumberOfClusters();
-    fTRDchi2=t->GetChi2();
-    for (Int_t i=0;i<fTRDncls;i++) index[i]=t->GetClusterIndex(i);
+    fTRDchi2  = t->GetChi2();
+    fTRDncls  = t->GetNumberOfClusters();
+    for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i);
+    
     fTRDsignal=t->GetPIDsignal();
     break;
   case kTRDbackup:
@@ -777,6 +1292,16 @@ AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
   return kTRUE;
 }
 
+void 
+AliESDtrack::SetOuterParam(const AliExternalTrackParam *p, ULong_t flags) {
+  //
+  // This is a direct setter for the outer track parameters
+  //
+  SetStatus(flags);
+  if (fOp) delete fOp;
+  fOp=new AliExternalTrackParam(*p);
+}
+
 Bool_t 
 AliESDtrack::GetOuterExternalParameters
                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
@@ -818,9 +1343,18 @@ Int_t AliESDtrack::GetNcls(Int_t idet) const
     ncls = fTRDncls;
     break;
   case 3:
-    if (fTOFindex != 0)
+    if (fTOFindex != -1)
       ncls = 1;
     break;
+  case 4: //PHOS
+    break;
+  case 5: //HMPID
+    if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
+      if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
+       ncls = 1;
+      }
+    }    
+    break;
   default:
     break;
   }
@@ -843,11 +1377,23 @@ Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
     ncls = GetTRDclusters(idx);
     break;
   case 3:
-    if (fTOFindex != 0) {
-      idx[0] = GetTOFcluster();
+    if (fTOFindex != -1) {
+      idx[0] = fTOFindex;
       ncls = 1;
     }
     break;
+  case 4: //PHOS
+    break;
+  case 5:
+    if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
+      if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
+       idx[0] = GetHMPIDcluIdx();
+       ncls = 1;
+      }
+    }    
+    break;
+  case 6: //EMCAL
+    break;
   default:
     break;
   }
@@ -886,11 +1432,71 @@ Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
   //---------------------------------------------------------------------
   if (idx!=0) {
      Int_t *index=fFriendTrack->GetITSindices();
-     for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) idx[i]=index[i];
+     for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) {
+         if ( (i>=fITSncls) && (i<6) ) idx[i]=-1;
+         else idx[i]=index[i];
+     }
   }
   return fITSncls;
 }
 
+//_______________________________________________________________________
+Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
+                                        Float_t &xloc,Float_t &zloc) const {
+  //----------------------------------------------------------------------
+  // This function encodes in the module number also the status of cluster association
+  // "status" can have the following values: 
+  // 1 "found" (cluster is associated), 
+  // 2 "dead" (module is dead from OCDB), 
+  // 3 "skipped" (module or layer forced to be skipped),
+  // 4 "outinz" (track out of z acceptance), 
+  // 5 "nocls" (no clusters in the road), 
+  // 6 "norefit" (cluster rejected during refit), 
+  // 7 "deadzspd" (holes in z in SPD)
+  // Also given are the coordinates of the crossing point of track and module
+  // (in the local module ref. system)
+  // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
+  //----------------------------------------------------------------------
+
+  if(fITSModule[ilayer]==-1) {
+    idet = -1;
+    status=0;
+    xloc=-99.; zloc=-99.;
+    return kFALSE;
+  }
+
+  Int_t module = fITSModule[ilayer];
+
+  idet = Int_t(module/1000000);
+
+  module -= idet*1000000;
+
+  status = Int_t(module/100000);
+
+  module -= status*100000;
+
+  Int_t signs = Int_t(module/10000);
+
+  module-=signs*10000;
+
+  Int_t xInt = Int_t(module/100);
+  module -= xInt*100;
+
+  Int_t zInt = module;
+
+  if(signs==1) { xInt*=1; zInt*=1; }
+  if(signs==2) { xInt*=1; zInt*=-1; }
+  if(signs==3) { xInt*=-1; zInt*=1; }
+  if(signs==4) { xInt*=-1; zInt*=-1; }
+
+  xloc = 0.1*(Float_t)xInt;
+  zloc = 0.1*(Float_t)zInt;
+
+  if(status==4) idet = -1;
+
+  return kTRUE;
+}
+
 //_______________________________________________________________________
 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
   //---------------------------------------------------------------------
@@ -947,6 +1553,18 @@ UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
   return fTRDncls;
 }
 
+//_______________________________________________________________________
+UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
+  //---------------------------------------------------------------------
+  // This function returns indices of the assigned TRD tracklets 
+  //---------------------------------------------------------------------
+  if (idx!=0) {
+     Int_t *index=fFriendTrack->GetTRDindices();
+     for (Int_t i=0; i<6/*AliESDfriendTrack::kMaxTRDcluster*/; i++) idx[i]=index[i];
+  }
+  return fTRDncls;
+}
+
 //_______________________________________________________________________
 void AliESDtrack::SetTRDpid(const Double_t *p) {  
   // Sets values for the probability of each particle type (in TRD)
@@ -973,6 +1591,61 @@ Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
   return fTRDr[iSpecies];
 }
 
+void  AliESDtrack::SetNumberOfTRDslices(Int_t n) {
+  //Sets the number of slices used for PID 
+  if (fTRDnSlices != 0) return;
+  fTRDnSlices=kTRDnPlanes*n;
+  fTRDslices=new Double32_t[fTRDnSlices];
+  for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=-1.;
+}
+
+void  AliESDtrack::SetTRDslice(Double_t q, Int_t plane, Int_t slice) {
+  //Sets the charge q in the slice of the plane
+  Int_t ns=GetNumberOfTRDslices();
+  if (ns==0) {
+    AliError("No TRD slices allocated for this track !");
+    return;
+  }
+
+  if ((plane<0) || (plane>=kTRDnPlanes)) {
+    AliError("Wrong TRD plane !");
+    return;
+  }
+  if ((slice<0) || (slice>=ns)) {
+    AliError("Wrong TRD slice !");
+    return;
+  }
+  Int_t n=plane*ns + slice;
+  fTRDslices[n]=q;
+}
+
+Double_t  AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const {
+  //Gets the charge from the slice of the plane
+  Int_t ns=GetNumberOfTRDslices();
+  if (ns==0) {
+    //AliError("No TRD slices allocated for this track !");
+    return -1.;
+  }
+
+  if ((plane<0) || (plane>=kTRDnPlanes)) {
+    AliError("Wrong TRD plane !");
+    return -1.;
+  }
+  if ((slice<-1) || (slice>=ns)) {
+    //AliError("Wrong TRD slice !");  
+    return -1.;
+  }
+
+  if (slice==-1) {
+    Double_t q=0.;
+    for (Int_t i=0; i<ns; i++) q+=fTRDslices[plane*ns + i];
+    return q/ns;
+  }
+
+  return fTRDslices[plane*ns + slice];
+}
+
+
 //_______________________________________________________________________
 void AliESDtrack::SetTOFpid(const Double_t *p) {  
   // Sets the probability of each particle type (in TOF)
@@ -1041,8 +1714,54 @@ void AliESDtrack::GetESDpid(Double_t *p) const {
 }
 
 //_______________________________________________________________________
-Bool_t AliESDtrack::RelateToVertex
-(const AliESDVertex *vtx, Double_t b, Double_t maxd) {
+Bool_t AliESDtrack::RelateToVertexTPC(const AliESDVertex *vtx, 
+Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
+  //
+  // Try to relate the TPC-only track parameters to the vertex "vtx", 
+  // if the (rough) transverse impact parameter is not bigger then "maxd". 
+  //            Magnetic field is "b" (kG).
+  //
+  // a) The TPC-only paramters are extapolated to the DCA to the vertex.
+  // b) The impact parameters and their covariance matrix are calculated.
+  // c) An attempt to constrain the TPC-only params to the vertex is done.
+  //    The constrained params are returned via "cParam".
+  //
+  // In the case of success, the returned value is kTRUE
+  // otherwise, it's kFALSE)
+  // 
+
+  if (!fTPCInner) return kFALSE;
+  if (!vtx) return kFALSE;
+
+  Double_t dz[2],cov[3];
+  if (!fTPCInner->PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
+
+  fdTPC = dz[0];
+  fzTPC = dz[1];  
+  fCddTPC = cov[0];
+  fCdzTPC = cov[1];
+  fCzzTPC = cov[2];
+  
+  Double_t covar[6]; vtx->GetCovMatrix(covar);
+  Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
+  Double_t c[3]={covar[2],0.,covar[5]};
+
+  Double_t chi2=GetPredictedChi2(p,c);
+  if (chi2>kVeryBig) return kFALSE;
+
+  fCchi2TPC=chi2;
+
+  if (!cParam) return kTRUE;
+
+  *cParam = *fTPCInner;
+  if (!cParam->Update(p,c)) return kFALSE;
+
+  return kTRUE;
+}
+
+//_______________________________________________________________________
+Bool_t AliESDtrack::RelateToVertex(const AliESDVertex *vtx, 
+Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
   //
   // Try to relate this track to the vertex "vtx", 
   // if the (rough) transverse impact parameter is not bigger then "maxd". 
@@ -1051,64 +1770,45 @@ Bool_t AliESDtrack::RelateToVertex
   // a) The track gets extapolated to the DCA to the vertex.
   // b) The impact parameters and their covariance matrix are calculated.
   // c) An attempt to constrain this track to the vertex is done.
+  //    The constrained params are returned via "cParam".
   //
-  //    In the case of success, the returned value is kTRUE
-  //    (otherwise, it's kFALSE)
+  // In the case of success, the returned value is kTRUE
+  // (otherwise, it's kFALSE)
   //  
 
   if (!vtx) return kFALSE;
 
-  Double_t alpha=GetAlpha();
-  Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
-  Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
-  Double_t xv= vtx->GetXv()*cs + vtx->GetYv()*sn;
-  Double_t yv=-vtx->GetXv()*sn + vtx->GetYv()*cs, zv=vtx->GetZv();
-  x-=xv; y-=yv;
-
-  //Estimate the impact parameter neglecting the track curvature
-  Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt(1.- snp*snp));
-  if (d > maxd) return kFALSE; 
-
-  //Propagate to the DCA
-  Double_t crv=kB2C*b*GetParameter()[4];
-  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
+  Double_t dz[2],cov[3];
+  if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
 
-  Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt(1.-snp*snp));
-  sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv);
-  if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
-  else cs=1.;
+  fD = dz[0];
+  fZ = dz[1];  
+  fCdd = cov[0];
+  fCdz = cov[1];
+  fCzz = cov[2];
+  
+  Double_t covar[6]; vtx->GetCovMatrix(covar);
+  Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
+  Double_t c[3]={covar[2],0.,covar[5]};
 
-  x = xv*cs + yv*sn;
-  yv=-xv*sn + yv*cs; xv=x;
+  Double_t chi2=GetPredictedChi2(p,c);
+  if (chi2>kVeryBig) return kFALSE;
 
-  if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
+  fCchi2=chi2;
 
-  fD = GetParameter()[0] - yv;
-  fZ = GetParameter()[1] - zv;
-  
-  Double_t cov[6]; vtx->GetCovMatrix(cov);
 
-  //***** Improvements by A.Dainese
-  alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
-  Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
-  fCdd = GetCovariance()[0] + s2ylocvtx;   // neglecting correlations
-  fCdz = GetCovariance()[1];               // between (x,y) and z
-  fCzz = GetCovariance()[2] + cov[5];      // in vertex's covariance matrix
-  //*****
+  //--- Could now these lines be removed ? ---
+  delete fCp;
+  fCp=new AliExternalTrackParam(*this);  
 
-  {//Try to constrain 
-    Double_t p[2]={yv,zv}, c[3]={cov[2],0.,cov[5]};
-    Double_t chi2=GetPredictedChi2(p,c);
+  if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
+  //----------------------------------------
 
-    if (chi2>77.) return kFALSE;
 
-    AliExternalTrackParam tmp(*this);
-    if (!tmp.Update(p,c)) return kFALSE;
+  if (!cParam) return kTRUE;
 
-    fCchi2=chi2;
-    if (!fCp) fCp=new AliExternalTrackParam();
-    new (fCp) AliExternalTrackParam(tmp);
-  }
+  *cParam = *this;
+  if (!cParam->Update(p,c)) return kFALSE; 
 
   return kTRUE;
 }
@@ -1116,7 +1816,7 @@ Bool_t AliESDtrack::RelateToVertex
 //_______________________________________________________________________
 void AliESDtrack::Print(Option_t *) const {
   // Prints info on the track
-  
+  AliExternalTrackParam::Print();
   printf("ESD track info\n") ; 
   Double_t p[AliPID::kSPECIESN] ; 
   Int_t index = 0 ; 
@@ -1139,7 +1839,7 @@ void AliESDtrack::Print(Option_t *) const {
     GetTRDpid(p) ; 
     for(index = 0 ; index < AliPID::kSPECIES; index++) 
       printf("%f, ", p[index]) ;
-    printf("\n           signal = %f\n", GetTRDsignal()) ;
+      printf("\n           signal = %f\n", GetTRDsignal()) ;
   }
   if( IsOn(kTOFpid) ){
     printf("From TOF: ") ; 
@@ -1157,3 +1857,59 @@ void AliESDtrack::Print(Option_t *) const {
   }
 } 
 
+
+//
+// Draw functionality
+// Origin: Marian Ivanov, Marian.Ivanov@cern.ch
+//
+void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
+  //
+  // Fill points in the polymarker
+  //
+  TObjArray arrayRef;
+  arrayRef.AddLast(new AliExternalTrackParam(*this));
+  if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
+  if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
+  //
+  Double_t mpos[3]={0,0,0};
+  Int_t entries=arrayRef.GetEntries();
+  for (Int_t i=0;i<entries;i++){
+    Double_t pos[3];
+    ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
+    mpos[0]+=pos[0]/entries;
+    mpos[1]+=pos[1]/entries;
+    mpos[2]+=pos[2]/entries;    
+  }
+  // Rotate to the mean position
+  //
+  Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
+  for (Int_t i=0;i<entries;i++){
+    Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
+    if (!res) delete arrayRef.RemoveAt(i);
+  }
+  Int_t counter=0;
+  for (Double_t r=minR; r<maxR; r+=stepR){
+    Double_t sweight=0;
+    Double_t mlpos[3]={0,0,0};
+    for (Int_t i=0;i<entries;i++){
+      Double_t point[3]={0,0,0};
+      AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
+      if (!param) continue;
+      if (param->GetXYZAt(r,magF,point)){
+       Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
+       sweight+=weight;
+       mlpos[0]+=point[0]*weight;
+       mlpos[1]+=point[1]*weight;
+       mlpos[2]+=point[2]*weight;
+      }
+    }
+    if (sweight>0){
+      mlpos[0]/=sweight;
+      mlpos[1]/=sweight;
+      mlpos[2]/=sweight;      
+      pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
+      printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
+      counter++;
+    }
+  }
+}