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=6;//t->GetNumberOfTracklets(); //t->GetNumberOfClusters();
 696     //for (Int_t i=0;i<fTRDncls;i++) index[i]=t->GetClusterIndex(i);
 697     for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i);
 698
 699     fTRDsignal=t->GetPIDsignal();
 700     break;
 701   case kTRDbackup:
 702     if (!fOp) fOp=new AliExternalTrackParam(*t);
 703     else
 704       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 705     fTRDncls0 = t->GetNumberOfClusters();
 706     break;
 707   case kTOFin:
 708     break;
 709   case kTOFout:
 710     break;
 711   case kTRDStop:
 712     break;
 713   default:
 714     AliError("Wrong flag !");
 715     return kFALSE;
 716   }
 717
 718   return rc;
 719 }
 720
 721 //_______________________________________________________________________
 722 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
 723   //---------------------------------------------------------------------
 724   // This function returns external representation of the track parameters
 725   //---------------------------------------------------------------------
 726   x=GetX();
 727   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
 728 }
 729
 730 //_______________________________________________________________________
 731 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
 732   //---------------------------------------------------------------------
 733   // This function returns external representation of the cov. matrix
 734   //---------------------------------------------------------------------
 735   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
 736 }
 737
 738 //_______________________________________________________________________
 739 Bool_t AliESDtrack::GetConstrainedExternalParameters
 740                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 741   //---------------------------------------------------------------------
 742   // This function returns the constrained external track parameters
 743   //---------------------------------------------------------------------
 744   if (!fCp) return kFALSE;
 745   alpha=fCp->GetAlpha();
 746   x=fCp->GetX();
 747   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
 748   return kTRUE;
 749 }
 750
 751 //_______________________________________________________________________
 752 Bool_t
 753 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
 754   //---------------------------------------------------------------------
 755   // This function returns the constrained external cov. matrix
 756   //---------------------------------------------------------------------
 757   if (!fCp) return kFALSE;
 758   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
 759   return kTRUE;
 760 }
 761
 762 Bool_t
 763 AliESDtrack::GetInnerExternalParameters
 764                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 765   //---------------------------------------------------------------------
 766   // This function returns external representation of the track parameters
 767   // at the inner layer of TPC
 768   //---------------------------------------------------------------------
 769   if (!fIp) return kFALSE;
 770   alpha=fIp->GetAlpha();
 771   x=fIp->GetX();
 772   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
 773   return kTRUE;
 774 }
 775
 776 Bool_t
 777 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
 778  //---------------------------------------------------------------------
 779  // This function returns external representation of the cov. matrix
 780  // at the inner layer of TPC
 781  //---------------------------------------------------------------------
 782   if (!fIp) return kFALSE;
 783   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
 784   return kTRUE;
 785 }
 786
 787 Bool_t
 788 AliESDtrack::GetOuterExternalParameters
 789                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 790   //---------------------------------------------------------------------
 791   // This function returns external representation of the track parameters
 792   // at the inner layer of TRD
 793   //---------------------------------------------------------------------
 794   if (!fOp) return kFALSE;
 795   alpha=fOp->GetAlpha();
 796   x=fOp->GetX();
 797   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
 798   return kTRUE;
 799 }
 800
 801 Bool_t
 802 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
 803  //---------------------------------------------------------------------
 804  // This function returns external representation of the cov. matrix
 805  // at the inner layer of TRD
 806  //---------------------------------------------------------------------
 807   if (!fOp) return kFALSE;
 808   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
 809   return kTRUE;
 810 }
 811
 812 Int_t AliESDtrack::GetNcls(Int_t idet) const
 813 {
 814   // Get number of clusters by subdetector index
 815   //
 816   Int_t ncls = 0;
 817   switch(idet){
 818   case 0:
 819     ncls = fITSncls;
 820     break;
 821   case 1:
 822     ncls = fTPCncls;
 823     break;
 824   case 2:
 825     ncls = fTRDncls;
 826     break;
 827   case 3:
 828     if (fTOFindex != 0)
 829       ncls = 1;
 830     break;
 831   default:
 832     break;
 833   }
 834   return ncls;
 835 }
 836
 837 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
 838 {
 839   // Get cluster index array by subdetector index
 840   //
 841   Int_t ncls = 0;
 842   switch(idet){
 843   case 0:
 844     ncls = GetITSclusters(idx);
 845     break;
 846   case 1:
 847     ncls = GetTPCclusters(idx);
 848     break;
 849   case 2:
 850     ncls = GetTRDclusters(idx);
 851     break;
 852   case 3:
 853     if (fTOFindex != 0) {
 854       idx[0] = GetTOFcluster();
 855       ncls = 1;
 856     }
 857     break;
 858   default:
 859     break;
 860   }
 861   return ncls;
 862 }
 863
 864 //_______________________________________________________________________
 865 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
 866   // Returns the array with integrated times for each particle hypothesis
 867   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
 868 }
 869
 870 //_______________________________________________________________________
 871 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
 872   // Sets the array with integrated times for each particle hypotesis
 873   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
 874 }
 875
 876 //_______________________________________________________________________
 877 void AliESDtrack::SetITSpid(const Double_t *p) {
 878   // Sets values for the probability of each particle type (in ITS)
 879   SetPIDValues(fITSr,p,AliPID::kSPECIES);
 880   SetStatus(AliESDtrack::kITSpid);
 881 }
 882
 883 //_______________________________________________________________________
 884 void AliESDtrack::GetITSpid(Double_t *p) const {
 885   // Gets the probability of each particle type (in ITS)
 886   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
 887 }
 888
 889 //_______________________________________________________________________
 890 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
 891   //---------------------------------------------------------------------
 892   // This function returns indices of the assgined ITS clusters
 893   //---------------------------------------------------------------------
 894   if (idx!=0) {
 895      Int_t *index=fFriendTrack->GetITSindices();
 896      for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) idx[i]=index[i];
 897   }
 898   return fITSncls;
 899 }
 900
 901 //_______________________________________________________________________
 902 Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
 903                                          Float_t &xloc,Float_t &zloc) const {
 904   //----------------------------------------------------------------------
 905   // This function encodes in the module number also the status of cluster association
 906   // "status" can have the following values:
 907   // 1 "found" (cluster is associated),
 908   // 2 "dead" (module is dead from OCDB),
 909   // 3 "skipped" (module or layer forced to be skipped),
 910   // 4 "outinz" (track out of z acceptance),
 911   // 5 "nocls" (no clusters in the road),
 912   // 6 "norefit" (cluster rejected during refit),
 913   // 7 "deadzspd" (holes in z in SPD)
 914   // Also given are the coordinates of the crossing point of track and module
 915   // (in the local module ref. system)
 916   // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
 917   //----------------------------------------------------------------------
 918
 919   if(fITSModule[ilayer]==-1) {
 920     AliError("fModule was not set !");
 921     idet = -1;
 922     status=0;
 923     xloc=-99.; zloc=-99.;
 924     return kFALSE;
 925   }
 926
 927   Int_t module = fITSModule[ilayer];
 928
 929   idet = Int_t(module/1000000);
 930
 931   module -= idet*1000000;
 932
 933   status = Int_t(module/100000);
 934
 935   module -= status*100000;
 936
 937   Int_t signs = Int_t(module/10000);
 938
 939   module-=signs*10000;
 940
 941   Int_t xInt = Int_t(module/100);
 942   module -= xInt*100;
 943
 944   Int_t zInt = module;
 945
 946   if(signs==1) { xInt*=1; zInt*=1; }
 947   if(signs==2) { xInt*=1; zInt*=-1; }
 948   if(signs==3) { xInt*=-1; zInt*=1; }
 949   if(signs==4) { xInt*=-1; zInt*=-1; }
 950
 951   xloc = 0.1*(Float_t)xInt;
 952   zloc = 0.1*(Float_t)zInt;
 953
 954   if(status==4) idet = -1;
 955
 956   return kTRUE;
 957 }
 958
 959 //_______________________________________________________________________
 960 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
 961   //---------------------------------------------------------------------
 962   // This function returns indices of the assgined ITS clusters
 963   //---------------------------------------------------------------------
 964   if (idx!=0) {
 965     Int_t *index=fFriendTrack->GetTPCindices();
 966     for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
 967   }
 968   return fTPCncls;
 969 }
 970
 971 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
 972   //
 973   // GetDensity of the clusters on given region between row0 and row1
 974   // Dead zone effect takin into acoount
 975   //
 976   Int_t good  = 0;
 977   Int_t found = 0;
 978   //
 979   Int_t *index=fFriendTrack->GetTPCindices();
 980   for (Int_t i=row0;i<=row1;i++){
 981     Int_t idx = index[i];
 982     if (idx!=-1)  good++;             // track outside of dead zone
 983     if (idx>0)    found++;
 984   }
 985   Float_t density=0.5;
 986   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
 987   return density;
 988 }
 989
 990 //_______________________________________________________________________
 991 void AliESDtrack::SetTPCpid(const Double_t *p) {
 992   // Sets values for the probability of each particle type (in TPC)
 993   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
 994   SetStatus(AliESDtrack::kTPCpid);
 995 }
 996
 997 //_______________________________________________________________________
 998 void AliESDtrack::GetTPCpid(Double_t *p) const {
 999   // Gets the probability of each particle type (in TPC)
1000   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
1001 }
1002
1003 //_______________________________________________________________________
1004 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
1005   //---------------------------------------------------------------------
1006   // This function returns indices of the assgined TRD clusters
1007   //---------------------------------------------------------------------
1008   if (idx!=0) {
1009      Int_t *index=fFriendTrack->GetTRDindices();
1010      for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
1011   }
1012   return fTRDncls;
1013 }
1014
1015 //_______________________________________________________________________
1016 UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
1017   //---------------------------------------------------------------------
1018   // This function returns indices of the assigned TRD tracklets
1019   //---------------------------------------------------------------------
1020   if (idx!=0) {
1021      Int_t *index=fFriendTrack->GetTRDindices();
1022      for (Int_t i=0; i<6/*AliESDfriendTrack::kMaxTRDcluster*/; i++) idx[i]=index[i];
1023   }
1024   return fTRDncls;
1025 }
1026
1027 //_______________________________________________________________________
1028 void AliESDtrack::SetTRDpid(const Double_t *p) {
1029   // Sets values for the probability of each particle type (in TRD)
1030   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
1031   SetStatus(AliESDtrack::kTRDpid);
1032 }
1033
1034 //_______________________________________________________________________
1035 void AliESDtrack::GetTRDpid(Double_t *p) const {
1036   // Gets the probability of each particle type (in TRD)
1037   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
1038 }
1039
1040 //_______________________________________________________________________
1041 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
1042 {
1043   // Sets the probability of particle type iSpecies to p (in TRD)
1044   fTRDr[iSpecies] = p;
1045 }
1046
1047 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
1048 {
1049   // Returns the probability of particle type iSpecies (in TRD)
1050   return fTRDr[iSpecies];
1051 }
1052
1053 //_______________________________________________________________________
1054 void AliESDtrack::SetTOFpid(const Double_t *p) {
1055   // Sets the probability of each particle type (in TOF)
1056   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
1057   SetStatus(AliESDtrack::kTOFpid);
1058 }
1059
1060 //_______________________________________________________________________
1061 void AliESDtrack::SetTOFLabel(const Int_t *p) {
1062   // Sets  (in TOF)
1063   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
1064 }
1065
1066 //_______________________________________________________________________
1067 void AliESDtrack::GetTOFpid(Double_t *p) const {
1068   // Gets probabilities of each particle type (in TOF)
1069   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
1070 }
1071
1072 //_______________________________________________________________________
1073 void AliESDtrack::GetTOFLabel(Int_t *p) const {
1074   // Gets (in TOF)
1075   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
1076 }
1077
1078 //_______________________________________________________________________
1079 void AliESDtrack::GetTOFInfo(Float_t *info) const {
1080   // Gets (in TOF)
1081   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
1082 }
1083
1084 //_______________________________________________________________________
1085 void AliESDtrack::SetTOFInfo(Float_t*info) {
1086   // Gets (in TOF)
1087   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
1088 }
1089
1090
1091
1092 //_______________________________________________________________________
1093 void AliESDtrack::SetHMPIDpid(const Double_t *p) {
1094   // Sets the probability of each particle type (in HMPID)
1095   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
1096   SetStatus(AliESDtrack::kHMPIDpid);
1097 }
1098
1099 //_______________________________________________________________________
1100 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
1101   // Gets probabilities of each particle type (in HMPID)
1102   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
1103 }
1104
1105
1106
1107 //_______________________________________________________________________
1108 void AliESDtrack::SetESDpid(const Double_t *p) {
1109   // Sets the probability of each particle type for the ESD track
1110   SetPIDValues(fR,p,AliPID::kSPECIES);
1111   SetStatus(AliESDtrack::kESDpid);
1112 }
1113
1114 //_______________________________________________________________________
1115 void AliESDtrack::GetESDpid(Double_t *p) const {
1116   // Gets probability of each particle type for the ESD track
1117   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1118 }
1119
1120 //_______________________________________________________________________
1121 Bool_t AliESDtrack::RelateToVertex
1122 (const AliESDVertex *vtx, Double_t b, Double_t maxd) {
1123   //
1124   // Try to relate this track to the vertex "vtx",
1125   // if the (rough) transverse impact parameter is not bigger then "maxd".
1126   //            Magnetic field is "b" (kG).
1127   //
1128   // a) The track gets extapolated to the DCA to the vertex.
1129   // b) The impact parameters and their covariance matrix are calculated.
1130   // c) An attempt to constrain this track to the vertex is done.
1131   //
1132   //    In the case of success, the returned value is kTRUE
1133   //    (otherwise, it's kFALSE)
1134   //
1135
1136   if (!vtx) return kFALSE;
1137
1138   Double_t dz[2],cov[3];
1139   if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
1140
1141   fD = dz[0];
1142   fZ = dz[1];
1143   fCdd = cov[0];
1144   fCdz = cov[1];
1145   fCzz = cov[2];
1146
1147   Double_t covar[6]; vtx->GetCovMatrix(covar);
1148   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
1149   Double_t c[3]={covar[2],0.,covar[5]};
1150
1151   Double_t chi2=GetPredictedChi2(p,c);
1152   if (chi2>77.) return kFALSE;
1153
1154   delete fCp;
1155   fCp=new AliExternalTrackParam(*this);
1156
1157   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
1158
1159   fCchi2=chi2;
1160   return kTRUE;
1161 }
1162
1163 //_______________________________________________________________________
1164 void AliESDtrack::Print(Option_t *) const {
1165   // Prints info on the track
1166
1167   printf("ESD track info\n") ;
1168   Double_t p[AliPID::kSPECIESN] ;
1169   Int_t index = 0 ;
1170   if( IsOn(kITSpid) ){
1171     printf("From ITS: ") ;
1172     GetITSpid(p) ;
1173     for(index = 0 ; index < AliPID::kSPECIES; index++)
1174       printf("%f, ", p[index]) ;
1175     printf("\n           signal = %f\n", GetITSsignal()) ;
1176   }
1177   if( IsOn(kTPCpid) ){
1178     printf("From TPC: ") ;
1179     GetTPCpid(p) ;
1180     for(index = 0 ; index < AliPID::kSPECIES; index++)
1181       printf("%f, ", p[index]) ;
1182     printf("\n           signal = %f\n", GetTPCsignal()) ;
1183   }
1184   if( IsOn(kTRDpid) ){
1185     printf("From TRD: ") ;
1186     GetTRDpid(p) ;
1187     for(index = 0 ; index < AliPID::kSPECIES; index++)
1188       printf("%f, ", p[index]) ;
1189     printf("\n           signal = %f\n", GetTRDsignal()) ;
1190   }
1191   if( IsOn(kTOFpid) ){
1192     printf("From TOF: ") ;
1193     GetTOFpid(p) ;
1194     for(index = 0 ; index < AliPID::kSPECIES; index++)
1195       printf("%f, ", p[index]) ;
1196     printf("\n           signal = %f\n", GetTOFsignal()) ;
1197   }
1198   if( IsOn(kHMPIDpid) ){
1199     printf("From HMPID: ") ;
1200     GetHMPIDpid(p) ;
1201     for(index = 0 ; index < AliPID::kSPECIES; index++)
1202       printf("%f, ", p[index]) ;
1203     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
1204   }
1205 }
1206
1207
1208 //
1209 // Draw functionality
1210 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
1211 //
1212 void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
1213   //
1214   // Fill points in the polymarker
1215   //
1216   TObjArray arrayRef;
1217   arrayRef.AddLast(new AliExternalTrackParam(*this));
1218   if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
1219   if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
1220   //
1221   Double_t mpos[3]={0,0,0};
1222   Int_t entries=arrayRef.GetEntries();
1223   for (Int_t i=0;i<entries;i++){
1224     Double_t pos[3];
1225     ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
1226     mpos[0]+=pos[0]/entries;
1227     mpos[1]+=pos[1]/entries;
1228     mpos[2]+=pos[2]/entries;
1229   }
1230   // Rotate to the mean position
1231   //
1232   Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
1233   for (Int_t i=0;i<entries;i++){
1234     Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
1235     if (!res) delete arrayRef.RemoveAt(i);
1236   }
1237   Int_t counter=0;
1238   for (Double_t r=minR; r<maxR; r+=stepR){
1239     Double_t sweight=0;
1240     Double_t mlpos[3]={0,0,0};
1241     for (Int_t i=0;i<entries;i++){
1242       Double_t point[3]={0,0,0};
1243       AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
1244       if (!param) continue;
1245       if (param->GetXYZAt(r,magF,point)){
1246         Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
1247         sweight+=weight;
1248         mlpos[0]+=point[0]*weight;
1249         mlpos[1]+=point[1]*weight;
1250         mlpos[2]+=point[2]*weight;
1251       }
1252     }
1253     if (sweight>0){
1254       mlpos[0]/=sweight;
1255       mlpos[1]/=sweight;
1256       mlpos[2]/=sweight;
1257       pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
1258       printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
1259       counter++;
1260     }
1261   }
1262 }