STEER/AliESDtrack.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 //-----------------------------------------------------------------
  16 //           Implementation of the ESD track class
  17 //   ESD = Event Summary Data
  18 //   This is the class to deal with during the phisics analysis of data
  19 //      Origin: Iouri Belikov, CERN
  20 //      e-mail: Jouri.Belikov@cern.ch
  21 //-----------------------------------------------------------------
  22
  23 #include <TMath.h>
  24 #include <TParticle.h>
  25
  26 #include "AliESDVertex.h"
  27 #include "AliESDtrack.h"
  28 #include "AliKalmanTrack.h"
  29 #include "AliLog.h"
  30 #include "AliTrackPointArray.h"
  31 #include "TPolyMarker3D.h"
  32
  33 ClassImp(AliESDtrack)
  34
  35 void SetPIDValues(Double_t * dest, const Double_t * src, Int_t n) {
  36   // This function copies "n" PID weights from "scr" to "dest"
  37   // and normalizes their sum to 1 thus producing conditional probabilities.
  38   // The negative weights are set to 0.
  39   // In case all the weights are non-positive they are replaced by
  40   // uniform probabilities
  41
  42   if (n<=0) return;
  43
  44   Float_t uniform = 1./(Float_t)n;
  45
  46   Float_t sum = 0;
  47   for (Int_t i=0; i<n; i++)
  48     if (src[i]>=0) {
  49       sum+=src[i];
  50       dest[i] = src[i];
  51     }
  52     else {
  53       dest[i] = 0;
  54     }
  55
  56   if(sum>0)
  57     for (Int_t i=0; i<n; i++) dest[i] /= sum;
  58   else
  59     for (Int_t i=0; i<n; i++) dest[i] = uniform;
  60 }
  61
  62 //_______________________________________________________________________
  63 AliESDtrack::AliESDtrack() :
  64   AliExternalTrackParam(),
  65   fCp(0),
  66   fIp(0),
  67   fTPCInner(0),
  68   fOp(0),
  69   fFriendTrack(new AliESDfriendTrack()),
  70   fTPCClusterMap(159),//number of padrows
  71   fTPCSharedMap(159),//number of padrows
  72   fFlags(0),
  73   fID(0),
  74   fLabel(0),
  75   fITSLabel(0),
  76   fTPCLabel(0),
  77   fTRDLabel(0),
  78   fTOFCalChannel(0),
  79   fTOFindex(0),
  80   fHMPIDqn(0),
  81   fHMPIDcluIdx(0),
  82   fEMCALindex(kEMCALNoMatch),
  83   fHMPIDtrkTheta(0),
  84   fHMPIDtrkPhi(0),
  85   fHMPIDsignal(0),
  86   fTrackLength(0),
  87   fD(0),fZ(0),
  88   fCdd(0),fCdz(0),fCzz(0),
  89   fCchi2(0),
  90   fITSchi2(0),
  91   fTPCchi2(0),
  92   fTRDchi2(0),
  93   fTOFchi2(0),
  94   fHMPIDchi2(0),
  95   fITSsignal(0),
  96   fTPCsignal(0),
  97   fTPCsignalS(0),
  98   fTRDsignal(0),
  99   fTRDQuality(0),
 100   fTRDBudget(0),
 101   fTOFsignal(0),
 102   fTOFsignalToT(0),
 103   fTOFsignalRaw(0),
 104   fTOFsignalDz(0),
 105   fHMPIDtrkX(0),
 106   fHMPIDtrkY(0),
 107   fHMPIDmipX(0),
 108   fHMPIDmipY(0),
 109   fTPCncls(0),
 110   fTPCnclsF(0),
 111   fTPCsignalN(0),
 112   fITSncls(0),
 113   fITSClusterMap(0),
 114   fTRDncls(0),
 115   fTRDncls0(0),
 116   fTRDpidQuality(0)
 117 {
 118   //
 119   // The default ESD constructor
 120   //
 121   Int_t i, j;
 122   for (i=0; i<AliPID::kSPECIES; i++) {
 123     fTrackTime[i]=0.;
 124     fR[i]=0.;
 125     fITSr[i]=0.;
 126     fTPCr[i]=0.;
 127     fTRDr[i]=0.;
 128     fTOFr[i]=0.;
 129     fHMPIDr[i]=0.;
 130   }
 131
 132   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 133   for (i=0; i<3; i++)   { fV0Indexes[i]=0;}
 134   for (i=0;i<kNPlane;i++) {
 135     for (j=0;j<kNSlice;j++) {
 136       fTRDsignals[i][j]=0.;
 137     }
 138     fTRDTimBin[i]=0;
 139   }
 140   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
 141   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
 142   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
 143   for (i=0;i<12;i++) {fITSModule[i]=-1;}
 144 }
 145
 146 //_______________________________________________________________________
 147 AliESDtrack::AliESDtrack(const AliESDtrack& track):
 148   AliExternalTrackParam(track),
 149   fCp(0),
 150   fIp(0),
 151   fTPCInner(0),
 152   fOp(0),
 153   fFriendTrack(0),
 154   fTPCClusterMap(track.fTPCClusterMap),
 155   fTPCSharedMap(track.fTPCSharedMap),
 156   fFlags(track.fFlags),
 157   fID(track.fID),
 158   fLabel(track.fLabel),
 159   fITSLabel(track.fITSLabel),
 160   fTPCLabel(track.fTPCLabel),
 161   fTRDLabel(track.fTRDLabel),
 162   fTOFCalChannel(track.fTOFCalChannel),
 163   fTOFindex(track.fTOFindex),
 164   fHMPIDqn(track.fHMPIDqn),
 165   fHMPIDcluIdx(track.fHMPIDcluIdx),
 166   fEMCALindex(track.fEMCALindex),
 167   fHMPIDtrkTheta(track.fHMPIDtrkTheta),
 168   fHMPIDtrkPhi(track.fHMPIDtrkPhi),
 169   fHMPIDsignal(track.fHMPIDsignal),
 170   fTrackLength(track.fTrackLength),
 171   fD(track.fD),fZ(track.fZ),
 172   fCdd(track.fCdd),fCdz(track.fCdz),fCzz(track.fCzz),
 173   fCchi2(track.fCchi2),
 174   fITSchi2(track.fITSchi2),
 175   fTPCchi2(track.fTPCchi2),
 176   fTRDchi2(track.fTRDchi2),
 177   fTOFchi2(track.fTOFchi2),
 178   fHMPIDchi2(track.fHMPIDchi2),
 179   fITSsignal(track.fITSsignal),
 180   fTPCsignal(track.fTPCsignal),
 181   fTPCsignalS(track.fTPCsignalS),
 182   fTRDsignal(track.fTRDsignal),
 183   fTRDQuality(track.fTRDQuality),
 184   fTRDBudget(track.fTRDBudget),
 185   fTOFsignal(track.fTOFsignal),
 186   fTOFsignalToT(track.fTOFsignalToT),
 187   fTOFsignalRaw(track.fTOFsignalRaw),
 188   fTOFsignalDz(track.fTOFsignalDz),
 189   fHMPIDtrkX(track.fHMPIDtrkX),
 190   fHMPIDtrkY(track.fHMPIDtrkY),
 191   fHMPIDmipX(track.fHMPIDmipX),
 192   fHMPIDmipY(track.fHMPIDmipY),
 193   fTPCncls(track.fTPCncls),
 194   fTPCnclsF(track.fTPCnclsF),
 195   fTPCsignalN(track.fTPCsignalN),
 196   fITSncls(track.fITSncls),
 197   fITSClusterMap(track.fITSClusterMap),
 198   fTRDncls(track.fTRDncls),
 199   fTRDncls0(track.fTRDncls0),
 200   fTRDpidQuality(track.fTRDpidQuality)
 201 {
 202   //
 203   //copy constructor
 204   //
 205   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i]=track.fTrackTime[i];
 206   for (Int_t i=0;i<AliPID::kSPECIES;i++)  fR[i]=track.fR[i];
 207   //
 208   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=track.fITSr[i];
 209   //
 210   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=track.fTPCr[i];
 211   for (Int_t i=0;i<4;i++) {fTPCPoints[i]=track.fTPCPoints[i];}
 212   for (Int_t i=0; i<3;i++)   { fKinkIndexes[i]=track.fKinkIndexes[i];}
 213   for (Int_t i=0; i<3;i++)   { fV0Indexes[i]=track.fV0Indexes[i];}
 214   //
 215   for (Int_t i=0;i<kNPlane;i++) {
 216     for (Int_t j=0;j<kNSlice;j++) {
 217       fTRDsignals[i][j]=track.fTRDsignals[i][j];
 218     }
 219     fTRDTimBin[i]=track.fTRDTimBin[i];
 220   }
 221   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i]=track.fTRDr[i];
 222   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i]=track.fTOFr[i];
 223   for (Int_t i=0;i<3;i++) fTOFLabel[i]=track.fTOFLabel[i];
 224   for (Int_t i=0;i<10;i++) fTOFInfo[i]=track.fTOFInfo[i];
 225   for (Int_t i=0;i<12;i++) fITSModule[i]=track.fITSModule[i];
 226   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i]=track.fHMPIDr[i];
 227
 228   if (track.fCp) fCp=new AliExternalTrackParam(*track.fCp);
 229   if (track.fIp) fIp=new AliExternalTrackParam(*track.fIp);
 230   if (track.fTPCInner) fTPCInner=new AliExternalTrackParam(*track.fTPCInner);
 231   if (track.fOp) fOp=new AliExternalTrackParam(*track.fOp);
 232
 233   if (track.fFriendTrack) fFriendTrack=new AliESDfriendTrack(*(track.fFriendTrack));
 234 }
 235
 236 //_______________________________________________________________________
 237 AliESDtrack::AliESDtrack(TParticle * part) :
 238   AliExternalTrackParam(),
 239   fCp(0),
 240   fIp(0),
 241   fTPCInner(0),
 242   fOp(0),
 243   fFriendTrack(0),
 244   fTPCClusterMap(159),//number of padrows
 245   fTPCSharedMap(159),//number of padrows
 246   fFlags(0),
 247   fID(0),
 248   fLabel(0),
 249   fITSLabel(0),
 250   fTPCLabel(0),
 251   fTRDLabel(0),
 252   fTOFCalChannel(0),
 253   fTOFindex(0),
 254   fHMPIDqn(0),
 255   fHMPIDcluIdx(0),
 256   fEMCALindex(kEMCALNoMatch),
 257   fHMPIDtrkTheta(0),
 258   fHMPIDtrkPhi(0),
 259   fHMPIDsignal(0),
 260   fTrackLength(0),
 261   fD(0),fZ(0),
 262   fCdd(0),fCdz(0),fCzz(0),
 263   fCchi2(0),
 264   fITSchi2(0),
 265   fTPCchi2(0),
 266   fTRDchi2(0),
 267   fTOFchi2(0),
 268   fHMPIDchi2(0),
 269   fITSsignal(0),
 270   fTPCsignal(0),
 271   fTPCsignalS(0),
 272   fTRDsignal(0),
 273   fTRDQuality(0),
 274   fTRDBudget(0),
 275   fTOFsignal(0),
 276   fTOFsignalToT(0),
 277   fTOFsignalRaw(0),
 278   fTOFsignalDz(0),
 279   fHMPIDtrkX(0),
 280   fHMPIDtrkY(0),
 281   fHMPIDmipX(0),
 282   fHMPIDmipY(0),
 283   fTPCncls(0),
 284   fTPCnclsF(0),
 285   fTPCsignalN(0),
 286   fITSncls(0),
 287   fITSClusterMap(0),
 288   fTRDncls(0),
 289   fTRDncls0(0),
 290   fTRDpidQuality(0)
 291 {
 292   //
 293   // ESD track from TParticle
 294   //
 295
 296   // Reset all the arrays
 297   Int_t i, j;
 298   for (i=0; i<AliPID::kSPECIES; i++) {
 299     fTrackTime[i]=0.;
 300     fR[i]=0.;
 301     fITSr[i]=0.;
 302     fTPCr[i]=0.;
 303     fTRDr[i]=0.;
 304     fTOFr[i]=0.;
 305     fHMPIDr[i]=0.;
 306   }
 307
 308   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 309   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
 310   for (i=0;i<kNPlane;i++) {
 311     for (j=0;j<kNSlice;j++) {
 312       fTRDsignals[i][j]=0.;
 313     }
 314     fTRDTimBin[i]=0;
 315   }
 316   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
 317   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
 318   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
 319   for (i=0;i<12;i++) {fITSModule[i]=-1;}
 320
 321   // Calculate the AliExternalTrackParam content
 322
 323   Double_t xref;
 324   Double_t alpha;
 325   Double_t param[5];
 326   Double_t covar[15];
 327
 328   // Calculate alpha: the rotation angle of the corresponding local system (TPC sector)
 329   alpha = part->Phi()*180./TMath::Pi();
 330   if (alpha<0) alpha+= 360.;
 331   if (alpha>360) alpha -= 360.;
 332
 333   Int_t sector = (Int_t)(alpha/20.);
 334   alpha = 10. + 20.*sector;
 335   alpha /= 180;
 336   alpha *= TMath::Pi();
 337
 338   // Covariance matrix: no errors, the parameters are exact
 339   for (i=0; i<15; i++) covar[i]=0.;
 340
 341   // Get the vertex of origin and the momentum
 342   TVector3 ver(part->Vx(),part->Vy(),part->Vz());
 343   TVector3 mom(part->Px(),part->Py(),part->Pz());
 344
 345   // Rotate to the local coordinate system (TPC sector)
 346   ver.RotateZ(-alpha);
 347   mom.RotateZ(-alpha);
 348
 349   // X of the referense plane
 350   xref = ver.X();
 351
 352   Int_t pdgCode = part->GetPdgCode();
 353
 354   Double_t charge =
 355     TDatabasePDG::Instance()->GetParticle(pdgCode)->Charge();
 356
 357   param[0] = ver.Y();
 358   param[1] = ver.Z();
 359   param[2] = TMath::Sin(mom.Phi());
 360   param[3] = mom.Pz()/mom.Pt();
 361   param[4] = TMath::Sign(1/mom.Pt(),charge);
 362
 363   // Set AliExternalTrackParam
 364   Set(xref, alpha, param, covar);
 365
 366   // Set the PID
 367   Int_t indexPID = 99;
 368
 369   switch (TMath::Abs(pdgCode)) {
 370
 371   case  11: // electron
 372     indexPID = 0;
 373     break;
 374
 375   case 13: // muon
 376     indexPID = 1;
 377     break;
 378
 379   case 211: // pion
 380     indexPID = 2;
 381     break;
 382
 383   case 321: // kaon
 384     indexPID = 3;
 385     break;
 386
 387   case 2212: // proton
 388     indexPID = 4;
 389     break;
 390
 391   default:
 392     break;
 393   }
 394
 395   // If the particle is not e,mu,pi,K or p the PID probabilities are set to 0
 396   if (indexPID < AliPID::kSPECIES) {
 397     fR[indexPID]=1.;
 398     fITSr[indexPID]=1.;
 399     fTPCr[indexPID]=1.;
 400     fTRDr[indexPID]=1.;
 401     fTOFr[indexPID]=1.;
 402     fHMPIDr[indexPID]=1.;
 403
 404   }
 405   // AliESD track label
 406   SetLabel(part->GetUniqueID());
 407
 408 }
 409
 410 //_______________________________________________________________________
 411 AliESDtrack::~AliESDtrack(){
 412   //
 413   // This is destructor according Coding Conventrions
 414   //
 415   //printf("Delete track\n");
 416   delete fIp;
 417   delete fTPCInner;
 418   delete fOp;
 419   delete fCp;
 420   delete fFriendTrack;
 421 }
 422
 423 void AliESDtrack::AddCalibObject(TObject * object){
 424   //
 425   // add calib object to the list
 426   //
 427   if (!fFriendTrack) fFriendTrack  = new AliESDfriendTrack;
 428   fFriendTrack->AddCalibObject(object);
 429 }
 430
 431 TObject *  AliESDtrack::GetCalibObject(Int_t index){
 432   //
 433   // return calib objct at given position
 434   //
 435   if (!fFriendTrack) return 0;
 436   return fFriendTrack->GetCalibObject(index);
 437 }
 438
 439
 440 //_______________________________________________________________________
 441 void AliESDtrack::MakeMiniESDtrack(){
 442   // Resets everything except
 443   // fFlags: Reconstruction status flags
 444   // fLabel: Track label
 445   // fID:  Unique ID of the track
 446   // fD: Impact parameter in XY-plane
 447   // fZ: Impact parameter in Z
 448   // fR[AliPID::kSPECIES]: combined "detector response probability"
 449   // Running track parameters in the base class (AliExternalTrackParam)
 450
 451   fTrackLength = 0;
 452
 453   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
 454
 455   // Reset track parameters constrained to the primary vertex
 456   delete fCp;fCp = 0;
 457   fCchi2 = 0;
 458
 459   // Reset track parameters at the inner wall of TPC
 460   delete fIp;fIp = 0;
 461   delete fTPCInner;fTPCInner=0;
 462   // Reset track parameters at the inner wall of the TRD
 463   delete fOp;fOp = 0;
 464
 465
 466   // Reset ITS track related information
 467   fITSchi2 = 0;
 468   fITSncls = 0;
 469   fITSClusterMap=0;
 470   fITSsignal = 0;
 471   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0;
 472   fITSLabel = 0;
 473
 474   // Reset TPC related track information
 475   fTPCchi2 = 0;
 476   fTPCncls = 0;
 477   fTPCnclsF = 0;
 478   fTPCClusterMap = 0;
 479   fTPCSharedMap = 0;
 480   fTPCsignal= 0;
 481   fTPCsignalS= 0;
 482   fTPCsignalN= 0;
 483   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0;
 484   fTPCLabel=0;
 485   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
 486   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
 487   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
 488
 489   // Reset TRD related track information
 490   fTRDchi2 = 0;
 491   fTRDncls = 0;
 492   fTRDncls0 = 0;
 493   fTRDsignal = 0;
 494   for (Int_t i=0;i<kNPlane;i++) {
 495     for (Int_t j=0;j<kNSlice;j++) {
 496       fTRDsignals[i][j] = 0;
 497     }
 498     fTRDTimBin[i]  = 0;
 499   }
 500   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0;
 501   fTRDLabel = 0;
 502   fTRDQuality  = 0;
 503   fTRDpidQuality = 0;
 504   fTRDBudget  = 0;
 505
 506   // Reset TOF related track information
 507   fTOFchi2 = 0;
 508   fTOFindex = 0;
 509   fTOFsignal = 0;
 510   fTOFCalChannel = 0;
 511   fTOFsignalToT = 0;
 512   fTOFsignalRaw = 0;
 513   fTOFsignalDz = 0;
 514   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
 515   for (Int_t i=0;i<3;i++) fTOFLabel[i] = 0;
 516   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
 517
 518   // Reset HMPID related track information
 519   fHMPIDchi2 = 0;
 520   fHMPIDqn = 0;
 521   fHMPIDcluIdx = 0;
 522   fHMPIDsignal = 0;
 523   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
 524   fHMPIDtrkTheta = 0;
 525   fHMPIDtrkPhi = 0;
 526   fHMPIDtrkX = 0;
 527   fHMPIDtrkY = 0;
 528   fHMPIDmipX = 0;
 529   fHMPIDmipY = 0;
 530   fEMCALindex = kEMCALNoMatch;
 531
 532   delete fFriendTrack; fFriendTrack = 0;
 533 }
 534 //_______________________________________________________________________
 535 Double_t AliESDtrack::GetMass() const {
 536   // Returns the mass of the most probable particle type
 537   Float_t max=0.;
 538   Int_t k=-1;
 539   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
 540     if (fR[i]>max) {k=i; max=fR[i];}
 541   }
 542   if (k==0) { // dE/dx "crossing points" in the TPC
 543      Double_t p=GetP();
 544      if ((p>0.38)&&(p<0.48))
 545         if (fR[0]<fR[3]*10.) return AliPID::ParticleMass(AliPID::kKaon);
 546      if ((p>0.75)&&(p<0.85))
 547         if (fR[0]<fR[4]*10.) return AliPID::ParticleMass(AliPID::kProton);
 548      return 0.00051;
 549   }
 550   if (k==1) return AliPID::ParticleMass(AliPID::kMuon);
 551   if (k==2||k==-1) return AliPID::ParticleMass(AliPID::kPion);
 552   if (k==3) return AliPID::ParticleMass(AliPID::kKaon);
 553   if (k==4) return AliPID::ParticleMass(AliPID::kProton);
 554   AliWarning("Undefined mass !");
 555   return AliPID::ParticleMass(AliPID::kPion);
 556 }
 557
 558 //______________________________________________________________________________
 559 Double_t AliESDtrack::E() const
 560 {
 561   // Returns the energy of the particle given its assumed mass.
 562   // Assumes the pion mass if the particle can't be identified properly.
 563
 564   Double_t m = M();
 565   Double_t p = P();
 566   return TMath::Sqrt(p*p + m*m);
 567 }
 568
 569 //______________________________________________________________________________
 570 Double_t AliESDtrack::Y() const
 571 {
 572   // Returns the rapidity of a particle given its assumed mass.
 573   // Assumes the pion mass if the particle can't be identified properly.
 574
 575   Double_t e = E();
 576   Double_t pz = Pz();
 577   if (e != TMath::Abs(pz)) { // energy was not equal to pz
 578     return 0.5*TMath::Log((e+pz)/(e-pz));
 579   } else { // energy was equal to pz
 580     return -999.;
 581   }
 582 }
 583
 584 //_______________________________________________________________________
 585 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
 586   //
 587   // This function updates track's running parameters
 588   //
 589   Int_t *index=0;
 590   Bool_t rc=kTRUE;
 591
 592   SetStatus(flags);
 593   fLabel=t->GetLabel();
 594
 595   if (t->IsStartedTimeIntegral()) {
 596     SetStatus(kTIME);
 597     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
 598     SetIntegratedLength(t->GetIntegratedLength());
 599   }
 600
 601   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 602
 603   switch (flags) {
 604
 605   case kITSin: case kITSout: case kITSrefit:
 606     fITSClusterMap=0;
 607     fITSncls=t->GetNumberOfClusters();
 608     index=fFriendTrack->GetITSindices();
 609     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
 610         index[i]=t->GetClusterIndex(i);
 611         if (i<fITSncls) {
 612            Int_t l=(index[i] & 0xf0000000) >> 28;
 613            SETBIT(fITSClusterMap,l);
 614         }
 615     }
 616     fITSchi2=t->GetChi2();
 617     fITSsignal=t->GetPIDsignal();
 618     fITSLabel = t->GetLabel();
 619     break;
 620
 621   case kTPCin: case kTPCrefit:
 622     fTPCLabel = t->GetLabel();
 623     if (flags==kTPCin)  fTPCInner=new AliExternalTrackParam(*t);
 624     if (!fIp) fIp=new AliExternalTrackParam(*t);
 625     else
 626       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 627   case kTPCout:
 628     index=fFriendTrack->GetTPCindices();
 629     if (flags & kTPCout){
 630       if (!fOp) fOp=new AliExternalTrackParam(*t);
 631       else
 632         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 633     }
 634     fTPCncls=t->GetNumberOfClusters();
 635     fTPCchi2=t->GetChi2();
 636
 637      {//prevrow must be declared in separate namespace, otherwise compiler cries:
 638       //"jump to case label crosses initialization of `Int_t prevrow'"
 639        Int_t prevrow = -1;
 640        //       for (Int_t i=0;i<fTPCncls;i++)
 641        for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++)
 642         {
 643           index[i]=t->GetClusterIndex(i);
 644           Int_t idx = index[i];
 645
 646           if (idx<0) continue;
 647
 648           // Piotr's Cluster Map for HBT
 649           // ### please change accordingly if cluster array is changing
 650           // to "New TPC Tracking" style (with gaps in array)
 651           Int_t sect = (idx&0xff000000)>>24;
 652           Int_t row = (idx&0x00ff0000)>>16;
 653           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
 654
 655           fTPCClusterMap.SetBitNumber(row,kTRUE);
 656
 657           //Fill the gap between previous row and this row with 0 bits
 658           //In case  ###  pleas change it as well - just set bit 0 in case there
 659           //is no associated clusters for current "i"
 660           if (prevrow < 0)
 661            {
 662              prevrow = row;//if previous bit was not assigned yet == this is the first one
 663            }
 664           else
 665            { //we don't know the order (inner to outer or reverse)
 666              //just to be save in case it is going to change
 667              Int_t n = 0, m = 0;
 668              if (prevrow < row)
 669               {
 670                 n = prevrow;
 671                 m = row;
 672               }
 673              else
 674               {
 675                 n = row;
 676                 m = prevrow;
 677               }
 678
 679              for (Int_t j = n+1; j < m; j++)
 680               {
 681                 fTPCClusterMap.SetBitNumber(j,kFALSE);
 682               }
 683              prevrow = row;
 684            }
 685           // End Of Piotr's Cluster Map for HBT
 686         }
 687      }
 688     fTPCsignal=t->GetPIDsignal();
 689     break;
 690
 691   case kTRDout: case kTRDin: case kTRDrefit:
 692     index     = fFriendTrack->GetTRDindices();
 693     fTRDLabel = t->GetLabel();
 694     fTRDchi2  = t->GetChi2();
 695     fTRDncls  = t->GetNumberOfClusters();
 696     for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i);
 697
 698     fTRDsignal=t->GetPIDsignal();
 699     break;
 700   case kTRDbackup:
 701     if (!fOp) fOp=new AliExternalTrackParam(*t);
 702     else
 703       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 704     fTRDncls0 = t->GetNumberOfClusters();
 705     break;
 706   case kTOFin:
 707     break;
 708   case kTOFout:
 709     break;
 710   case kTRDStop:
 711     break;
 712   default:
 713     AliError("Wrong flag !");
 714     return kFALSE;
 715   }
 716
 717   return rc;
 718 }
 719
 720 //_______________________________________________________________________
 721 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
 722   //---------------------------------------------------------------------
 723   // This function returns external representation of the track parameters
 724   //---------------------------------------------------------------------
 725   x=GetX();
 726   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
 727 }
 728
 729 //_______________________________________________________________________
 730 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
 731   //---------------------------------------------------------------------
 732   // This function returns external representation of the cov. matrix
 733   //---------------------------------------------------------------------
 734   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
 735 }
 736
 737 //_______________________________________________________________________
 738 Bool_t AliESDtrack::GetConstrainedExternalParameters
 739                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 740   //---------------------------------------------------------------------
 741   // This function returns the constrained external track parameters
 742   //---------------------------------------------------------------------
 743   if (!fCp) return kFALSE;
 744   alpha=fCp->GetAlpha();
 745   x=fCp->GetX();
 746   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
 747   return kTRUE;
 748 }
 749
 750 //_______________________________________________________________________
 751 Bool_t
 752 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
 753   //---------------------------------------------------------------------
 754   // This function returns the constrained external cov. matrix
 755   //---------------------------------------------------------------------
 756   if (!fCp) return kFALSE;
 757   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
 758   return kTRUE;
 759 }
 760
 761 Bool_t
 762 AliESDtrack::GetInnerExternalParameters
 763                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 764   //---------------------------------------------------------------------
 765   // This function returns external representation of the track parameters
 766   // at the inner layer of TPC
 767   //---------------------------------------------------------------------
 768   if (!fIp) return kFALSE;
 769   alpha=fIp->GetAlpha();
 770   x=fIp->GetX();
 771   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
 772   return kTRUE;
 773 }
 774
 775 Bool_t
 776 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
 777  //---------------------------------------------------------------------
 778  // This function returns external representation of the cov. matrix
 779  // at the inner layer of TPC
 780  //---------------------------------------------------------------------
 781   if (!fIp) return kFALSE;
 782   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
 783   return kTRUE;
 784 }
 785
 786 Bool_t
 787 AliESDtrack::GetOuterExternalParameters
 788                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 789   //---------------------------------------------------------------------
 790   // This function returns external representation of the track parameters
 791   // at the inner layer of TRD
 792   //---------------------------------------------------------------------
 793   if (!fOp) return kFALSE;
 794   alpha=fOp->GetAlpha();
 795   x=fOp->GetX();
 796   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
 797   return kTRUE;
 798 }
 799
 800 Bool_t
 801 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
 802  //---------------------------------------------------------------------
 803  // This function returns external representation of the cov. matrix
 804  // at the inner layer of TRD
 805  //---------------------------------------------------------------------
 806   if (!fOp) return kFALSE;
 807   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
 808   return kTRUE;
 809 }
 810
 811 Int_t AliESDtrack::GetNcls(Int_t idet) const
 812 {
 813   // Get number of clusters by subdetector index
 814   //
 815   Int_t ncls = 0;
 816   switch(idet){
 817   case 0:
 818     ncls = fITSncls;
 819     break;
 820   case 1:
 821     ncls = fTPCncls;
 822     break;
 823   case 2:
 824     ncls = fTRDncls;
 825     break;
 826   case 3:
 827     if (fTOFindex != 0)
 828       ncls = 1;
 829     break;
 830   default:
 831     break;
 832   }
 833   return ncls;
 834 }
 835
 836 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
 837 {
 838   // Get cluster index array by subdetector index
 839   //
 840   Int_t ncls = 0;
 841   switch(idet){
 842   case 0:
 843     ncls = GetITSclusters(idx);
 844     break;
 845   case 1:
 846     ncls = GetTPCclusters(idx);
 847     break;
 848   case 2:
 849     ncls = GetTRDclusters(idx);
 850     break;
 851   case 3:
 852     if (fTOFindex != 0) {
 853       idx[0] = GetTOFcluster();
 854       ncls = 1;
 855     }
 856     break;
 857   default:
 858     break;
 859   }
 860   return ncls;
 861 }
 862
 863 //_______________________________________________________________________
 864 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
 865   // Returns the array with integrated times for each particle hypothesis
 866   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
 867 }
 868
 869 //_______________________________________________________________________
 870 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
 871   // Sets the array with integrated times for each particle hypotesis
 872   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
 873 }
 874
 875 //_______________________________________________________________________
 876 void AliESDtrack::SetITSpid(const Double_t *p) {
 877   // Sets values for the probability of each particle type (in ITS)
 878   SetPIDValues(fITSr,p,AliPID::kSPECIES);
 879   SetStatus(AliESDtrack::kITSpid);
 880 }
 881
 882 //_______________________________________________________________________
 883 void AliESDtrack::GetITSpid(Double_t *p) const {
 884   // Gets the probability of each particle type (in ITS)
 885   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
 886 }
 887
 888 //_______________________________________________________________________
 889 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
 890   //---------------------------------------------------------------------
 891   // This function returns indices of the assgined ITS clusters
 892   //---------------------------------------------------------------------
 893   if (idx!=0) {
 894      Int_t *index=fFriendTrack->GetITSindices();
 895      for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) idx[i]=index[i];
 896   }
 897   return fITSncls;
 898 }
 899
 900 //_______________________________________________________________________
 901 Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
 902                                          Float_t &xloc,Float_t &zloc) const {
 903   //----------------------------------------------------------------------
 904   // This function encodes in the module number also the status of cluster association
 905   // "status" can have the following values:
 906   // 1 "found" (cluster is associated),
 907   // 2 "dead" (module is dead from OCDB),
 908   // 3 "skipped" (module or layer forced to be skipped),
 909   // 4 "outinz" (track out of z acceptance),
 910   // 5 "nocls" (no clusters in the road),
 911   // 6 "norefit" (cluster rejected during refit),
 912   // 7 "deadzspd" (holes in z in SPD)
 913   // Also given are the coordinates of the crossing point of track and module
 914   // (in the local module ref. system)
 915   // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
 916   //----------------------------------------------------------------------
 917
 918   if(fITSModule[ilayer]==-1) {
 919     AliError("fModule was not set !");
 920     idet = -1;
 921     status=0;
 922     xloc=-99.; zloc=-99.;
 923     return kFALSE;
 924   }
 925
 926   Int_t module = fITSModule[ilayer];
 927
 928   idet = Int_t(module/1000000);
 929
 930   module -= idet*1000000;
 931
 932   status = Int_t(module/100000);
 933
 934   module -= status*100000;
 935
 936   Int_t signs = Int_t(module/10000);
 937
 938   module-=signs*10000;
 939
 940   Int_t xInt = Int_t(module/100);
 941   module -= xInt*100;
 942
 943   Int_t zInt = module;
 944
 945   if(signs==1) { xInt*=1; zInt*=1; }
 946   if(signs==2) { xInt*=1; zInt*=-1; }
 947   if(signs==3) { xInt*=-1; zInt*=1; }
 948   if(signs==4) { xInt*=-1; zInt*=-1; }
 949
 950   xloc = 0.1*(Float_t)xInt;
 951   zloc = 0.1*(Float_t)zInt;
 952
 953   if(status==4) idet = -1;
 954
 955   return kTRUE;
 956 }
 957
 958 //_______________________________________________________________________
 959 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
 960   //---------------------------------------------------------------------
 961   // This function returns indices of the assgined ITS clusters
 962   //---------------------------------------------------------------------
 963   if (idx!=0) {
 964     Int_t *index=fFriendTrack->GetTPCindices();
 965     for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
 966   }
 967   return fTPCncls;
 968 }
 969
 970 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
 971   //
 972   // GetDensity of the clusters on given region between row0 and row1
 973   // Dead zone effect takin into acoount
 974   //
 975   Int_t good  = 0;
 976   Int_t found = 0;
 977   //
 978   Int_t *index=fFriendTrack->GetTPCindices();
 979   for (Int_t i=row0;i<=row1;i++){
 980     Int_t idx = index[i];
 981     if (idx!=-1)  good++;             // track outside of dead zone
 982     if (idx>0)    found++;
 983   }
 984   Float_t density=0.5;
 985   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
 986   return density;
 987 }
 988
 989 //_______________________________________________________________________
 990 void AliESDtrack::SetTPCpid(const Double_t *p) {
 991   // Sets values for the probability of each particle type (in TPC)
 992   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
 993   SetStatus(AliESDtrack::kTPCpid);
 994 }
 995
 996 //_______________________________________________________________________
 997 void AliESDtrack::GetTPCpid(Double_t *p) const {
 998   // Gets the probability of each particle type (in TPC)
 999   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
1000 }
1001
1002 //_______________________________________________________________________
1003 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
1004   //---------------------------------------------------------------------
1005   // This function returns indices of the assgined TRD clusters
1006   //---------------------------------------------------------------------
1007   if (idx!=0) {
1008      Int_t *index=fFriendTrack->GetTRDindices();
1009      for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
1010   }
1011   return fTRDncls;
1012 }
1013
1014 //_______________________________________________________________________
1015 UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
1016   //---------------------------------------------------------------------
1017   // This function returns indices of the assigned TRD tracklets
1018   //---------------------------------------------------------------------
1019   if (idx!=0) {
1020      Int_t *index=fFriendTrack->GetTRDindices();
1021      for (Int_t i=0; i<6/*AliESDfriendTrack::kMaxTRDcluster*/; i++) idx[i]=index[i];
1022   }
1023   return fTRDncls;
1024 }
1025
1026 //_______________________________________________________________________
1027 void AliESDtrack::SetTRDpid(const Double_t *p) {
1028   // Sets values for the probability of each particle type (in TRD)
1029   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
1030   SetStatus(AliESDtrack::kTRDpid);
1031 }
1032
1033 //_______________________________________________________________________
1034 void AliESDtrack::GetTRDpid(Double_t *p) const {
1035   // Gets the probability of each particle type (in TRD)
1036   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
1037 }
1038
1039 //_______________________________________________________________________
1040 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
1041 {
1042   // Sets the probability of particle type iSpecies to p (in TRD)
1043   fTRDr[iSpecies] = p;
1044 }
1045
1046 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
1047 {
1048   // Returns the probability of particle type iSpecies (in TRD)
1049   return fTRDr[iSpecies];
1050 }
1051
1052 //_______________________________________________________________________
1053 void AliESDtrack::SetTOFpid(const Double_t *p) {
1054   // Sets the probability of each particle type (in TOF)
1055   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
1056   SetStatus(AliESDtrack::kTOFpid);
1057 }
1058
1059 //_______________________________________________________________________
1060 void AliESDtrack::SetTOFLabel(const Int_t *p) {
1061   // Sets  (in TOF)
1062   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
1063 }
1064
1065 //_______________________________________________________________________
1066 void AliESDtrack::GetTOFpid(Double_t *p) const {
1067   // Gets probabilities of each particle type (in TOF)
1068   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
1069 }
1070
1071 //_______________________________________________________________________
1072 void AliESDtrack::GetTOFLabel(Int_t *p) const {
1073   // Gets (in TOF)
1074   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
1075 }
1076
1077 //_______________________________________________________________________
1078 void AliESDtrack::GetTOFInfo(Float_t *info) const {
1079   // Gets (in TOF)
1080   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
1081 }
1082
1083 //_______________________________________________________________________
1084 void AliESDtrack::SetTOFInfo(Float_t*info) {
1085   // Gets (in TOF)
1086   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
1087 }
1088
1089
1090
1091 //_______________________________________________________________________
1092 void AliESDtrack::SetHMPIDpid(const Double_t *p) {
1093   // Sets the probability of each particle type (in HMPID)
1094   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
1095   SetStatus(AliESDtrack::kHMPIDpid);
1096 }
1097
1098 //_______________________________________________________________________
1099 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
1100   // Gets probabilities of each particle type (in HMPID)
1101   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
1102 }
1103
1104
1105
1106 //_______________________________________________________________________
1107 void AliESDtrack::SetESDpid(const Double_t *p) {
1108   // Sets the probability of each particle type for the ESD track
1109   SetPIDValues(fR,p,AliPID::kSPECIES);
1110   SetStatus(AliESDtrack::kESDpid);
1111 }
1112
1113 //_______________________________________________________________________
1114 void AliESDtrack::GetESDpid(Double_t *p) const {
1115   // Gets probability of each particle type for the ESD track
1116   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1117 }
1118
1119 //_______________________________________________________________________
1120 Bool_t AliESDtrack::RelateToVertex
1121 (const AliESDVertex *vtx, Double_t b, Double_t maxd) {
1122   //
1123   // Try to relate this track to the vertex "vtx",
1124   // if the (rough) transverse impact parameter is not bigger then "maxd".
1125   //            Magnetic field is "b" (kG).
1126   //
1127   // a) The track gets extapolated to the DCA to the vertex.
1128   // b) The impact parameters and their covariance matrix are calculated.
1129   // c) An attempt to constrain this track to the vertex is done.
1130   //
1131   //    In the case of success, the returned value is kTRUE
1132   //    (otherwise, it's kFALSE)
1133   //
1134
1135   if (!vtx) return kFALSE;
1136
1137   Double_t dz[2],cov[3];
1138   if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
1139
1140   fD = dz[0];
1141   fZ = dz[1];
1142   fCdd = cov[0];
1143   fCdz = cov[1];
1144   fCzz = cov[2];
1145
1146   Double_t covar[6]; vtx->GetCovMatrix(covar);
1147   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
1148   Double_t c[3]={covar[2],0.,covar[5]};
1149
1150   Double_t chi2=GetPredictedChi2(p,c);
1151   if (chi2>77.) return kFALSE;
1152
1153   delete fCp;
1154   fCp=new AliExternalTrackParam(*this);
1155
1156   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
1157
1158   fCchi2=chi2;
1159   return kTRUE;
1160 }
1161
1162 //_______________________________________________________________________
1163 void AliESDtrack::Print(Option_t *) const {
1164   // Prints info on the track
1165
1166   printf("ESD track info\n") ;
1167   Double_t p[AliPID::kSPECIESN] ;
1168   Int_t index = 0 ;
1169   if( IsOn(kITSpid) ){
1170     printf("From ITS: ") ;
1171     GetITSpid(p) ;
1172     for(index = 0 ; index < AliPID::kSPECIES; index++)
1173       printf("%f, ", p[index]) ;
1174     printf("\n           signal = %f\n", GetITSsignal()) ;
1175   }
1176   if( IsOn(kTPCpid) ){
1177     printf("From TPC: ") ;
1178     GetTPCpid(p) ;
1179     for(index = 0 ; index < AliPID::kSPECIES; index++)
1180       printf("%f, ", p[index]) ;
1181     printf("\n           signal = %f\n", GetTPCsignal()) ;
1182   }
1183   if( IsOn(kTRDpid) ){
1184     printf("From TRD: ") ;
1185     GetTRDpid(p) ;
1186     for(index = 0 ; index < AliPID::kSPECIES; index++)
1187       printf("%f, ", p[index]) ;
1188     printf("\n           signal = %f\n", GetTRDsignal()) ;
1189   }
1190   if( IsOn(kTOFpid) ){
1191     printf("From TOF: ") ;
1192     GetTOFpid(p) ;
1193     for(index = 0 ; index < AliPID::kSPECIES; index++)
1194       printf("%f, ", p[index]) ;
1195     printf("\n           signal = %f\n", GetTOFsignal()) ;
1196   }
1197   if( IsOn(kHMPIDpid) ){
1198     printf("From HMPID: ") ;
1199     GetHMPIDpid(p) ;
1200     for(index = 0 ; index < AliPID::kSPECIES; index++)
1201       printf("%f, ", p[index]) ;
1202     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
1203   }
1204 }
1205
1206
1207 //
1208 // Draw functionality
1209 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
1210 //
1211 void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
1212   //
1213   // Fill points in the polymarker
1214   //
1215   TObjArray arrayRef;
1216   arrayRef.AddLast(new AliExternalTrackParam(*this));
1217   if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
1218   if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
1219   //
1220   Double_t mpos[3]={0,0,0};
1221   Int_t entries=arrayRef.GetEntries();
1222   for (Int_t i=0;i<entries;i++){
1223     Double_t pos[3];
1224     ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
1225     mpos[0]+=pos[0]/entries;
1226     mpos[1]+=pos[1]/entries;
1227     mpos[2]+=pos[2]/entries;
1228   }
1229   // Rotate to the mean position
1230   //
1231   Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
1232   for (Int_t i=0;i<entries;i++){
1233     Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
1234     if (!res) delete arrayRef.RemoveAt(i);
1235   }
1236   Int_t counter=0;
1237   for (Double_t r=minR; r<maxR; r+=stepR){
1238     Double_t sweight=0;
1239     Double_t mlpos[3]={0,0,0};
1240     for (Int_t i=0;i<entries;i++){
1241       Double_t point[3]={0,0,0};
1242       AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
1243       if (!param) continue;
1244       if (param->GetXYZAt(r,magF,point)){
1245         Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
1246         sweight+=weight;
1247         mlpos[0]+=point[0]*weight;
1248         mlpos[1]+=point[1]*weight;
1249         mlpos[2]+=point[2]*weight;
1250       }
1251     }
1252     if (sweight>0){
1253       mlpos[0]/=sweight;
1254       mlpos[1]/=sweight;
1255       mlpos[2]/=sweight;
1256       pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
1257       printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
1258       counter++;
1259     }
1260   }
1261 }