[u/mrichter/AliRoot.git] / ANALYSIS / AliTrackPoints.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

////////////////////////////////////////////////////////////
//                                                        //
// class AliTrackPoints                                //
//                                                        //
// used by Anti-Merging cut                               //
// contains set of poits the lay on track trajectory      //
// according to reconstructed track parameters -          //
// NOT CLUSTERS POSITIONS!!!                              //
// Anti-Merging cut is applied only on tracks coming from //
// different events (that are use to fill deniminators)   //
//                                                        //
////////////////////////////////////////////////////////////

#include "AliTrackPoints.h"
#include "AliLog.h"

#include <TClonesArray.h>
#include <TFile.h>
#include <TMath.h>

#include "AliESDtrack.h"
#include "AliTPCtrack.h"
#include "AliTrackReference.h"
#include "AliITStrackV2.h"

#include "AliRun.h"
#include "AliMagF.h"
#include "AliTracker.h"
#include "AliESD.h"
#include "AliRunLoader.h"
#include "AliTPCtrack.h"
#include "TTree.h"
#include "TBranch.h"
#include "TH2D.h"

ClassImp(AliTrackPoints)

Int_t AliTrackPoints::fgDebug = 0;

AliTrackPoints::AliTrackPoints():
 fN(0),
 fX(0x0),
 fY(0x0),
 fZ(0x0)
{
  //constructor
}

/***************************************************************/

AliTrackPoints::AliTrackPoints(AliTrackPoints::ETypes type, AliESDtrack* track, Float_t mf):
 fN(0),
 fX(0x0),
 fY(0x0),
 fZ(0x0)
{
  //constructor 
  //tupe -  what kind of track points should be calculated
  //mf - magnetic field in [kG] = [T]*10.0
  switch (type)
   {
     case kITS:
       fN = 6;
       fX = new Float_t[fN];
       fY = new Float_t[fN];
       fZ = new Float_t[fN];
       MakeITSPoints(track);
       break;

     case kITSInnerFromVertexOuterFromTPC:
       fN = 6;
       fX = new Float_t[fN];
       fY = new Float_t[fN];
       fZ = new Float_t[fN];
       MakeITSPointsInnerFromVertexOuterFromTPC(track,mf);
       break;

     default:
       Info("AliTrackPoints","Not recognized type");
   }
   
}

/***************************************************************/

AliTrackPoints::AliTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr, Float_t r0):
 fN(n),
 fX(new Float_t[fN]),
 fY(new Float_t[fN]),
 fZ(new Float_t[fN])
{
  //constructor
  //mf - magnetic field in kG - needed to calculate curvature out of Pt
  //r0 - starting radius
  //dr - calculate points every dr cm, default every 30cm
  if (track == 0x0)
   {
     Error("AliTrackPoints","ESD track is null");
     fN = 0;
     delete [] fX;
     delete [] fY;
     delete [] fZ;
     fX = fY = fZ = 0x0;
     return;
   }

  if ( ((track->GetStatus() & AliESDtrack::kTPCrefit) == kFALSE)&& 
       ((track->GetStatus() & AliESDtrack::kTPCin) == kFALSE)  )
   {
     //could happend: its stand alone tracking
     AliDebug(3,"This ESD track does not contain TPC information");
       
     fN = 0;
     delete [] fX;
     delete [] fY;
     delete [] fZ;
     fX = fY = fZ = 0x0;
     
     return;
   }
  
  Double_t alpha,x;
  Double_t par[5];
  track->GetInnerExternalParameters(alpha,x,par);//get properties of the track
  if (par[4] == 0)
   {
     Error("AliTrackPoints","This ESD track seem not to contain TPC information (curv is 0)");
     return;
   }

  if (mf == 0.0)
   {
     Error("AliTrackPoints","Zero Magnetic field passed as parameter.");
     return;
   }
   
  //Double_t alpha = track->GetInnerAlpha();
  Double_t cc = 1000./0.299792458/mf;//conversion constant
  Double_t c=par[4]/cc;
  
  MakePoints(dr,r0,x,par,c,alpha);
  
}

/***************************************************************/

AliTrackPoints::AliTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr, Float_t r0):
 fN(n),
 fX(new Float_t[fN]),
 fY(new Float_t[fN]),
 fZ(new Float_t[fN])
{
  //constructor
  //r0 starting radius
  //dr - calculate points every dr cm, default every 30cm
  if (track == 0x0)
   {
     Error("AliTrackPoints","TPC track is null");
     fN = 0;
     delete [] fX;
     delete [] fY;
     delete [] fZ;
     fX = fY = fZ = 0x0;
     return;
   }
  track->PropagateTo(r0);

 //* This formation is now fixed in the following way:                         *
 //*      external param0:   local Y-coordinate of a track (cm)                *
 //*      external param1:   local Z-coordinate of a track (cm)                *
 //*      external param2:   local sine of the track momentum azimuth angle    *
 //*      external param3:   tangent of the track momentum dip angle           *
 //*      external param4:   1/pt (1/(GeV/c))                                  *
    
  Double_t x = 0;
  Double_t par[5];
  track->GetExternalParameters(x,par);     //get properties of the track
  
  Double_t alpha = track->GetAlpha();
  Double_t c=track->GetC();
  MakePoints(dr,r0,x,par,c,alpha);
}

/***************************************************************/

AliTrackPoints::~AliTrackPoints()
{
  //destructor
  delete [] fX;
  delete [] fY;
  delete [] fZ;
}

/***************************************************************/

void AliTrackPoints::MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha)
{
  //Calculates points starting at radius r0
  //spacing every dr (in radial direction)
  // according to track parameters
  // x - position in sector local reference frame. x i parallel to R and sector is symmetric with respect to x
  // par - track external parameters; array with 5 elements;  look at AliTPCtrack.h or AliESDtrack.h for their meaning
  // c - track curvature
  // alpha - sector's rotation angle (phi) == angle needed for local to global transformation
  
  Double_t y = par[0];
  Double_t z0 = par[1];
      
  Double_t phi0local = TMath::ATan2(y,x);
  Double_t phi0global = phi0local + alpha;
  
  if (phi0local<0) phi0local+=2*TMath::Pi();
  if (phi0local>=2.*TMath::Pi()) phi0local-=2*TMath::Pi();
  
  if (phi0global<0) phi0global+=2*TMath::Pi();
  if (phi0global>=2.*TMath::Pi()) phi0global-=2*TMath::Pi();
  
  Double_t r = TMath::Hypot(x,y);
  

  AliDebug(9,Form("Radius0 %f, Real Radius %f",r0,r)); 
  
  AliDebug(5,Form("Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local));
    
  Double_t eta = x*c - par[2] ;//par[2] = fX*C - eta; eta==fP2 ; C==fP4
  
  //this calculattions are assuming new (current) model 
  Double_t tmp = par[2];
  tmp = 1. - tmp*tmp;
  tmp = c*y + TMath::Sqrt(tmp);
  Double_t dca=(TMath::Hypot(eta,tmp) - 1. )/TMath::Abs(c);

  //Here is old model Cold=Cnew/2.
  Double_t dcasq = dca*dca;
  Double_t c2 = c/2.;
  Double_t cst1 = (1.+c2*dca)*dca;//first constant
  Double_t cst2 = 1. + 2.*c2*dca;//second constant
      
  Double_t factorPhi0 = TMath::ASin((c2*r + cst1/r)/cst2);
  Double_t factorZ0 = TMath::ASin(c2*TMath::Sqrt((r*r-dcasq)/cst2))*par[3]/c2;
 
  for(Int_t i = 0; i<fN; i++)
   {
     Double_t rc = r0 + i*dr;
     Double_t ftmp = (c2*rc + cst1/rc)/cst2;
     if (ftmp > 1.0) 
      {
        AliDebug(1,Form("ASin argument > 1 %f:",ftmp));
        ftmp=1.0;
      }
     else if (ftmp < -1.0) 
      {
        AliDebug(1,Form("ASin argument < -1 %f:",ftmp));
        ftmp=-1.0;
      }
      
     Double_t factorPhi = TMath::ASin( ftmp );//factor phi od rc
     Double_t phi = phi0global + factorPhi - factorPhi0;
     
     ftmp = (rc*rc-dcasq)/cst2;
     if (ftmp < 0.0) 
      {
        AliDebug(1,Form("Sqrt argument < 0: %f",ftmp));
        ftmp=0.0;
      }
     
     ftmp = c2*TMath::Sqrt(ftmp);
     if (ftmp > 1.0) 
      {
        AliDebug(1,Form("ASin argument > 1: %f",ftmp));
        ftmp=1.0;
      }
     else if (ftmp < -1.0) 
      {
        AliDebug(2,Form("ASin argument < -1: %f",ftmp));
        ftmp=-1.0;
      }
     Double_t factorZ = TMath::ASin(ftmp)*par[3]/c2;
     fZ[i] = z0 + factorZ - factorZ0;
     fX[i] = rc*TMath::Cos(phi);
     fY[i] = rc*TMath::Sin(phi);
     
     AliDebug(3,Form("AliTrackPoints","X %f Y %f Z %f R asked %f R obtained %f",
             fX[i],fY[i],fZ[i],rc,TMath::Hypot(fX[i],fY[i])));
   }
}

/***************************************************************/

void AliTrackPoints::MakeITSPoints(AliESDtrack* track)
{
//Calculates points in ITS
// z=R*Pz/Pt
 AliITStrackV2 itstrack(*track,kTRUE);
 Double_t x,y,z;
 static const Double_t kR[6] = {4.0, 7.0, 14.9, 23.8, 39.1, 43.6};
 for (Int_t i = 0; i < 6; i++)
  {
    itstrack.GetGlobalXYZat(kR[i],x,y,z);
    fX[i] = x;
    fY[i] = y;
    fZ[i] = z;
//    Info("MakeITSPoints","X %f Y %f Z %f R asked %f R obtained %f",
//             fX[i],fY[i],fZ[i],r[i],TMath::Hypot(fX[i],fY[i]));
  }   
 
}

/***************************************************************/

void AliTrackPoints::MakeITSPointsInnerFromVertexOuterFromTPC(AliESDtrack* track, Float_t mf)
{
//makes trackpoints for ITS  
//for 3 inner layers calculates out of the vector at vertex
//for 3 outer ---------------//------------------ at inner TPC  

 static const Double_t kR[6] = {4.0, 7.0, 14.9, 23.8, 39.1, 43.6};
 AliITStrackV2 itstrack(*track,kTRUE);
 Double_t x,y,z;
 for (Int_t i = 0; i < 3; i++)
  {
    itstrack.GetGlobalXYZat(kR[i],x,y,z);
    fX[i] = x;
    fY[i] = y;
    fZ[i] = z;
    AliDebug(3,Form("X %f Y %f Z %f R asked %f R obtained %f",
                fX[i],fY[i],fZ[i],kR[i],TMath::Hypot(fX[i],fY[i])));
  }   
 
 for (Int_t i = 3; i < 6; i++)
  {
    Float_t ax,ay,az;
    AliTrackPoints tmptp(1,track,mf,0,kR[i]);
    tmptp.PositionAt(0,ax,ay,az);
    fX[i] = ax;
    fY[i] = ay;
    fZ[i] = az;
    AliDebug(3,Form("X %f Y %f Z %f R asked %f R obtained %f",
                fX[i],fY[i],fZ[i],kR[i],TMath::Hypot(fX[i],fY[i])));
  }
 
}
    
/***************************************************************/

void AliTrackPoints::PositionAt(Int_t n, Float_t &x,Float_t &y,Float_t &z)
{
  //returns position at point n
  if ((n<0) || (n>=fN))
   {
     Error("PositionAt","Point %d out of range",n);
     return;
   }
  
  x = fX[n];
  y = fY[n];
  z = fZ[n];
  AliDebug(2,Form("n %d; X %f; Y %f; Z %f",n,x,y,z));

}

/***************************************************************/

void AliTrackPoints::Move(Float_t x, Float_t y, Float_t z)
{
//Moves all points about vector
 for (Int_t i = 0; i<fN; i++)
   {
     fX[i]+=x;
     fY[i]+=y;
     fZ[i]+=z;
   }   
}

/***************************************************************/

Double_t AliTrackPoints::AvarageDistance(const AliTrackPoints& tr)
{
  //returns the aritmethic avarage distance between two tracks
//  Info("AvarageDistance","Entered");
  if ( (fN <= 0) || (tr.fN <=0) )
   {
     AliDebug(1,"One of tracks is empty");
     return -1;
   }

  if (fN != tr.fN)
   {
     Warning("AvarageDistance","Number of points is not equal");
     return -1;
   }
   
  Double_t sum = 0;
  for (Int_t i = 0; i<fN; i++)
   {
     AliDebug(10,Form("radii: %f %f",TMath::Hypot(fX[i],fY[i]),TMath::Hypot(tr.fX[i],tr.fY[i])));
//       Float_t r1sq = fX[i]*fX[i]+fY[i]*fY[i];
//       Float_t r2sq = tr.fX[i]*tr.fX[i]+tr.fY[i]*tr.fY[i];

      
     Double_t dx = fX[i]-tr.fX[i];
     Double_t dy = fY[i]-tr.fY[i];
     Double_t dz = fZ[i]-tr.fZ[i];
     sum+=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
     
     AliDebug(2,Form("Diff: x ,y z: %f , %f, %f",dx,dy,dz));
     AliDebug(2,Form("xxyyzz %f %f %f %f %f %f",
            fX[i],tr.fX[i],fY[i],tr.fY[i],fZ[i],tr.fZ[i]));
   }
   
  Double_t retval = sum/((Double_t)fN);
  AliDebug(1,Form("Avarage distance is %f.",retval));

  return retval;
}

/***************************************************************/

void AliTrackPoints::Print(Option_t* /*option*/) const
{
  // Prints the coordinates of the track points
  Info("Print","There is %d points",fN);
  for(Int_t i = 0; i < fN; i++)
   {
     Info("Print","%d: %f %f %f",i,fX[i],fY[i],fZ[i]);
   }

}

/***************************************************************/

void AliTrackPoints::Testesd(Int_t entr,const char* fname )
{
  // This is a former macro converted to member function.
  // It is used to test the functionality of the class using ESD tracks
  delete gAlice;
  gAlice = 0x0;
  AliRunLoader* rl = AliRunLoader::Open();
  rl->LoadgAlice();
  
  Float_t mf = rl->GetAliRun()->Field()->SolenoidField();
 
  
  TFile* fFile = TFile::Open(fname);
  
  if (fFile == 0x0)
   {
     printf("testesd: There is no suche a ESD file\n");
     return;
   }
  AliESD* esd = dynamic_cast<AliESD*>(fFile->Get("0"));
  AliESDtrack *t = esd->GetTrack(entr);
  if (t == 0x0)
   {
     ::Error("testesd","Can not get track %d",entr);
     return;
   }

  
  Int_t nTrackPoints = 170;
  AliTrackPoints* tp = new AliTrackPoints(nTrackPoints,t,mf,1.);
  
  Float_t xmin = -250;
  Float_t xmax =  250;
  
  Float_t ymin = -250;
  Float_t ymax =  250;

  Float_t zmin = -250;
  Float_t zmax =  250;
  
  TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
  TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
  TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);

  TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
  TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
  TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);

  hxyt->SetDirectory(0x0);   
  hxy->SetDirectory(0x0);   
  hxyTR->SetDirectory(0x0);   
  
  hxzt->SetDirectory(0x0);   
  hxz->SetDirectory(0x0);   
  hxzTR->SetDirectory(0x0);   

  Float_t x,y,z;
   
  for (Int_t i = 0;i<nTrackPoints;i++)
   {  
      Double_t r = 84.1+i;
      tp->PositionAt(i,x,y,z);
      hxy->Fill(x,y);
      hxz->Fill(x,z);
      printf("Rdemanded %f\n",r);
      printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
      
   }

  rl->LoadTrackRefs();
  TTree* treeTR = rl->TreeTR();
  TBranch* b = treeTR->GetBranch("TPC");
  
  TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
  AliTrackReference* tref;
  b->SetAddress(&trackrefs);
  
  Int_t tlab = TMath::Abs(t->GetLabel());

  Int_t netr = (Int_t)treeTR->GetEntries();
  printf("Found %d entries in TR tree\n",netr);
  
  for (Int_t e = 0; e < netr; e++)
   { 
     treeTR->GetEntry(e);
     tref = (AliTrackReference*)trackrefs->At(0);
     if (tref == 0x0) continue;
     if (tref->GetTrack() != tlab) continue;
    
     printf("Found %d entries in TR array\n",trackrefs->GetEntries());
     
     for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
      {
        tref = (AliTrackReference*)trackrefs->At(i);
        if (tref->GetTrack() != tlab) continue;
        x = tref->X();
        y = tref->Y();
        z = tref->Z();
        printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());

        hxzTR->Fill(x,z);
        hxyTR->Fill(x,y);
        for (Int_t j = 1; j < 10; j++)
         {
            hxyTR->Fill(x, y+j*0.1);
            hxyTR->Fill(x, y-j*0.1);
            hxyTR->Fill(x+j*0.1,y);
            hxyTR->Fill(x-j*0.1,y);

            hxzTR->Fill(x,z-j*0.1);
            hxzTR->Fill(x,z+j*0.1);
            hxzTR->Fill(x-j*0.1,z);
            hxzTR->Fill(x+j*0.1,z);
         }
      }
      break; 
    }  
  hxy->Draw("");
//  hxzt->Draw("same");
  hxyTR->Draw("same");
  
  delete rl;
}

/***************************************************************/

void AliTrackPoints::Testtpc(Int_t entr)
{
  // This is a former macro converted to member function.
  // It is used to test the functionality of the class using TPC tracks
  delete gAlice;
  gAlice = 0x0;
  AliRunLoader* rl = AliRunLoader::Open();
  AliLoader* l = rl->GetLoader("TPCLoader");
  rl->LoadgAlice();
  AliTracker::SetFieldMap(rl->GetAliRun()->Field(),kTRUE);
  l->LoadTracks();
  AliTPCtrack* t = new AliTPCtrack();
  TBranch* b=l->TreeT()->GetBranch("tracks");
  b->SetAddress(&t);
  l->TreeT()->GetEntry(entr);  
  Int_t nTrackPoints = 160;
  AliTrackPoints* tp = new AliTrackPoints(nTrackPoints,t,1.);
  
  Float_t xmin = -250;
  Float_t xmax =  250;
  
  Float_t ymin = -250;
  Float_t ymax =  250;

  Float_t zmin = -250;
  Float_t zmax =  250;
  
  TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
  TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
  TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);

  TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
  TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
  TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);

  hxyt->SetDirectory(0x0);   
  hxy->SetDirectory(0x0);   
  hxyTR->SetDirectory(0x0);   
  
  hxzt->SetDirectory(0x0);   
  hxz->SetDirectory(0x0);   
  hxzTR->SetDirectory(0x0);   

  Float_t x,y,z;
   
  for (Int_t i = 0;i<nTrackPoints;i++)
   {  
      Double_t r = 84.1+i;
      tp->PositionAt(i,x,y,z);
      hxy->Fill(x,y);
      hxz->Fill(x,z);
      printf("Rdemanded %f\n",r);
      printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));

      //BUT they are local!!!!
      t->PropagateTo(r);
//      Double_t phi = t->Phi();
      Double_t rl = TMath::Hypot(t->GetX(),t->GetY());//real radius 
      
      Double_t alpha = t->GetAlpha();
      Double_t salpha = TMath::Sin(alpha);
      Double_t calpha = TMath::Cos(alpha);
      x = t->GetX()*calpha - t->GetY()*salpha;
      y = t->GetX()*salpha + t->GetY()*calpha;
      z = t->GetZ();
      
      printf("tx %f ty %f tz %f Rt = %f R from XY %f\n",x,y,z,TMath::Hypot(x,y),rl);
      
      printf("tpz - tz %f\n",z-t->GetZ());
      printf("\n");
      hxyt->Fill(x,y);
      hxzt->Fill(x,z);
      
   }

  rl->LoadTrackRefs();
  TTree* treeTR = rl->TreeTR();
  b = treeTR->GetBranch("TPC");
  
  TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
  AliTrackReference* tref;
  b->SetAddress(&trackrefs);
  
  Int_t tlab = TMath::Abs(t->GetLabel());

  Int_t netr = (Int_t)treeTR->GetEntries();
  printf("Found %d entries in TR tree\n",netr);
  
  for (Int_t e = 0; e < netr; e++)
   { 
     treeTR->GetEntry(e);
     tref = (AliTrackReference*)trackrefs->At(0);
     if (tref == 0x0) continue;
     if (tref->GetTrack() != tlab) continue;
    
     printf("Found %d entries in TR array\n",trackrefs->GetEntries());
     
     for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
      {
        tref = (AliTrackReference*)trackrefs->At(i);
        if (tref->GetTrack() != tlab) continue;
        x = tref->X();
        y = tref->Y();
        z = tref->Z();
        printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());

        hxzTR->Fill(x,z);
        hxyTR->Fill(x,y);
        for (Int_t j = 1; j < 10; j++)
         {
            hxyTR->Fill(x, y+j*0.1);
            hxyTR->Fill(x, y-j*0.1);
            hxyTR->Fill(x+j*0.1,y);
            hxyTR->Fill(x-j*0.1,y);

            hxzTR->Fill(x,z-j*0.1);
            hxzTR->Fill(x,z+j*0.1);
            hxzTR->Fill(x-j*0.1,z);
            hxzTR->Fill(x+j*0.1,z);
         }
      }
      break; 
    }  
  hxz->Draw("");
//  hxzt->Draw("same");
  hxzTR->Draw("same");
  
  delete rl;
}
Commit	Line	Data
c7ffd78f	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	/* $Id$ */
	17
a5556ea5	18	////////////////////////////////////////////////////////////
	19	// //
	20	// class AliTrackPoints //
	21	// //
	22	// used by Anti-Merging cut //
	23	// contains set of poits the lay on track trajectory //
	24	// according to reconstructed track parameters - //
	25	// NOT CLUSTERS POSITIONS!!! //
	26	// Anti-Merging cut is applied only on tracks coming from //
	27	// different events (that are use to fill deniminators) //
	28	// //
	29	////////////////////////////////////////////////////////////
	30
0206ddfb	31	#include "AliTrackPoints.h"
	32	#include "AliLog.h"
	33
a5556ea5	34	#include <TClonesArray.h>
	35	#include <TFile.h>
	36	#include <TMath.h>
	37
	38	#include "AliESDtrack.h"
a5556ea5	39	#include "AliTPCtrack.h"
a5556ea5	40	#include "AliTrackReference.h"
78d7c6d3	41	#include "AliITStrackV2.h"
a5556ea5	42
0206ddfb	43	#include "AliRun.h"
7d0f8548	44	#include "AliMagF.h"
7d0f8548	45	#include "AliTracker.h"
0206ddfb	46	#include "AliESD.h"
	47	#include "AliRunLoader.h"
	48	#include "AliTPCtrack.h"
	49	#include "TTree.h"
	50	#include "TBranch.h"
	51	#include "TH2D.h"
	52
a5556ea5	53	ClassImp(AliTrackPoints)
	54
	55	Int_t AliTrackPoints::fgDebug = 0;
78d7c6d3	56
a5556ea5	57	AliTrackPoints::AliTrackPoints():
	58	fN(0),
	59	fX(0x0),
	60	fY(0x0),
	61	fZ(0x0)
	62	{
	63	//constructor
	64	}
0206ddfb	65
a5556ea5	66	/***************************************************************/
a5556ea5	67
0f0355ab	68	AliTrackPoints::AliTrackPoints(AliTrackPoints::ETypes type, AliESDtrack* track, Float_t mf):
78d7c6d3	69	fN(0),
	70	fX(0x0),
	71	fY(0x0),
	72	fZ(0x0)
	73	{
0f0355ab	74	//constructor
	75	//tupe - what kind of track points should be calculated
	76	//mf - magnetic field in [kG] = [T]*10.0
78d7c6d3	77	switch (type)
	78	{
	79	case kITS:
78d7c6d3	80	fN = 6;
	81	fX = new Float_t[fN];
	82	fY = new Float_t[fN];
	83	fZ = new Float_t[fN];
	84	MakeITSPoints(track);
	85	break;
	86
0f0355ab	87	case kITSInnerFromVertexOuterFromTPC:
	88	fN = 6;
	89	fX = new Float_t[fN];
	90	fY = new Float_t[fN];
	91	fZ = new Float_t[fN];
	92	MakeITSPointsInnerFromVertexOuterFromTPC(track,mf);
	93	break;
	94
78d7c6d3	95	default:
	96	Info("AliTrackPoints","Not recognized type");
	97	}
	98
	99	}
0206ddfb	100
78d7c6d3	101	/***************************************************************/
78d7c6d3	102
a5556ea5	103	AliTrackPoints::AliTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr, Float_t r0):
	104	fN(n),
	105	fX(new Float_t[fN]),
	106	fY(new Float_t[fN]),
	107	fZ(new Float_t[fN])
	108	{
	109	//constructor
0f0355ab	110	//mf - magnetic field in kG - needed to calculate curvature out of Pt
a5556ea5	111	//r0 - starting radius
	112	//dr - calculate points every dr cm, default every 30cm
	113	if (track == 0x0)
	114	{
	115	Error("AliTrackPoints","ESD track is null");
	116	fN = 0;
	117	delete [] fX;
	118	delete [] fY;
	119	delete [] fZ;
	120	fX = fY = fZ = 0x0;
	121	return;
	122	}
	123
	124	if ( ((track->GetStatus() & AliESDtrack::kTPCrefit) == kFALSE)&&
	125	((track->GetStatus() & AliESDtrack::kTPCin) == kFALSE) )
	126	{
	127	//could happend: its stand alone tracking
c7ffd78f	128	AliDebug(3,"This ESD track does not contain TPC information");
a5556ea5	129
	130	fN = 0;
	131	delete [] fX;
	132	delete [] fY;
	133	delete [] fZ;
	134	fX = fY = fZ = 0x0;
	135
	136	return;
	137	}
	138
c0b978f0	139	Double_t alpha,x;
a5556ea5	140	Double_t par[5];
c0b978f0	141	track->GetInnerExternalParameters(alpha,x,par);//get properties of the track
a5556ea5	142	if (par[4] == 0)
	143	{
	144	Error("AliTrackPoints","This ESD track seem not to contain TPC information (curv is 0)");
	145	return;
	146	}
	147
	148	if (mf == 0.0)
	149	{
	150	Error("AliTrackPoints","Zero Magnetic field passed as parameter.");
	151	return;
	152	}
	153
c0b978f0	154	//Double_t alpha = track->GetInnerAlpha();
a5556ea5	155	Double_t cc = 1000./0.299792458/mf;//conversion constant
	156	Double_t c=par[4]/cc;
	157
	158	MakePoints(dr,r0,x,par,c,alpha);
	159
	160	}
0206ddfb	161
a5556ea5	162	/***************************************************************/
	163
	164	AliTrackPoints::AliTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr, Float_t r0):
	165	fN(n),
	166	fX(new Float_t[fN]),
	167	fY(new Float_t[fN]),
	168	fZ(new Float_t[fN])
	169	{
	170	//constructor
	171	//r0 starting radius
	172	//dr - calculate points every dr cm, default every 30cm
	173	if (track == 0x0)
	174	{
	175	Error("AliTrackPoints","TPC track is null");
	176	fN = 0;
	177	delete [] fX;
	178	delete [] fY;
	179	delete [] fZ;
	180	fX = fY = fZ = 0x0;
	181	return;
	182	}
	183	track->PropagateTo(r0);
	184
	185	//* This formation is now fixed in the following way: *
	186	//* external param0: local Y-coordinate of a track (cm) *
	187	//* external param1: local Z-coordinate of a track (cm) *
	188	//* external param2: local sine of the track momentum azimuth angle *
	189	//* external param3: tangent of the track momentum dip angle *
	190	//* external param4: 1/pt (1/(GeV/c)) *
	191
	192	Double_t x = 0;
	193	Double_t par[5];
	194	track->GetExternalParameters(x,par); //get properties of the track
	195
	196	Double_t alpha = track->GetAlpha();
	197	Double_t c=track->GetC();
	198	MakePoints(dr,r0,x,par,c,alpha);
0206ddfb	199	}
0206ddfb	200
78d7c6d3	201	/***************************************************************/
	202
	203	AliTrackPoints::~AliTrackPoints()
	204	{
	205	//destructor
	206	delete [] fX;
	207	delete [] fY;
	208	delete [] fZ;
	209	}
0206ddfb	210
78d7c6d3	211	/***************************************************************/
a5556ea5	212
	213	void AliTrackPoints::MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha)
	214	{
	215	//Calculates points starting at radius r0
	216	//spacing every dr (in radial direction)
	217	// according to track parameters
	218	// x - position in sector local reference frame. x i parallel to R and sector is symmetric with respect to x
	219	// par - track external parameters; array with 5 elements; look at AliTPCtrack.h or AliESDtrack.h for their meaning
	220	// c - track curvature
	221	// alpha - sector's rotation angle (phi) == angle needed for local to global transformation
	222
	223	Double_t y = par[0];
	224	Double_t z0 = par[1];
	225
	226	Double_t phi0local = TMath::ATan2(y,x);
	227	Double_t phi0global = phi0local + alpha;
	228
	229	if (phi0local<0) phi0local+=2*TMath::Pi();
	230	if (phi0local>=2.TMath::Pi()) phi0local-=2TMath::Pi();
	231
	232	if (phi0global<0) phi0global+=2*TMath::Pi();
	233	if (phi0global>=2.TMath::Pi()) phi0global-=2TMath::Pi();
	234
	235	Double_t r = TMath::Hypot(x,y);
	236
	237
c7ffd78f	238	AliDebug(9,Form("Radius0 %f, Real Radius %f",r0,r));
a5556ea5	239
c7ffd78f	240	AliDebug(5,Form("Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local));
a5556ea5	241
	242	Double_t eta = xc - par[2] ;//par[2] = fXC - eta; eta==fP2 ; C==fP4
	243
	244	//this calculattions are assuming new (current) model
	245	Double_t tmp = par[2];
	246	tmp = 1. - tmp*tmp;
	247	tmp = c*y + TMath::Sqrt(tmp);
	248	Double_t dca=(TMath::Hypot(eta,tmp) - 1. )/TMath::Abs(c);
	249
	250	//Here is old model Cold=Cnew/2.
	251	Double_t dcasq = dca*dca;
	252	Double_t c2 = c/2.;
	253	Double_t cst1 = (1.+c2dca)dca;//first constant
	254	Double_t cst2 = 1. + 2.c2dca;//second constant
	255
	256	Double_t factorPhi0 = TMath::ASin((c2*r + cst1/r)/cst2);
	257	Double_t factorZ0 = TMath::ASin(c2TMath::Sqrt((rr-dcasq)/cst2))*par[3]/c2;
	258
	259	for(Int_t i = 0; i<fN; i++)
	260	{
	261	Double_t rc = r0 + i*dr;
	262	Double_t ftmp = (c2*rc + cst1/rc)/cst2;
	263	if (ftmp > 1.0)
	264	{
c7ffd78f	265	AliDebug(1,Form("ASin argument > 1 %f:",ftmp));
a5556ea5	266	ftmp=1.0;
	267	}
	268	else if (ftmp < -1.0)
	269	{
c7ffd78f	270	AliDebug(1,Form("ASin argument < -1 %f:",ftmp));
a5556ea5	271	ftmp=-1.0;
	272	}
	273
	274	Double_t factorPhi = TMath::ASin( ftmp );//factor phi od rc
	275	Double_t phi = phi0global + factorPhi - factorPhi0;
	276
	277	ftmp = (rc*rc-dcasq)/cst2;
	278	if (ftmp < 0.0)
	279	{
c7ffd78f	280	AliDebug(1,Form("Sqrt argument < 0: %f",ftmp));
a5556ea5	281	ftmp=0.0;
	282	}
	283
	284	ftmp = c2*TMath::Sqrt(ftmp);
	285	if (ftmp > 1.0)
	286	{
c7ffd78f	287	AliDebug(1,Form("ASin argument > 1: %f",ftmp));
a5556ea5	288	ftmp=1.0;
	289	}
	290	else if (ftmp < -1.0)
	291	{
c7ffd78f	292	AliDebug(2,Form("ASin argument < -1: %f",ftmp));
a5556ea5	293	ftmp=-1.0;
	294	}
	295	Double_t factorZ = TMath::ASin(ftmp)*par[3]/c2;
	296	fZ[i] = z0 + factorZ - factorZ0;
	297	fX[i] = rc*TMath::Cos(phi);
	298	fY[i] = rc*TMath::Sin(phi);
	299
c7ffd78f	300	AliDebug(3,Form("AliTrackPoints","X %f Y %f Z %f R asked %f R obtained %f",
c7ffd78f	301	fX[i],fY[i],fZ[i],rc,TMath::Hypot(fX[i],fY[i])));
a5556ea5	302	}
a5556ea5	303	}
0206ddfb	304
a5556ea5	305	/***************************************************************/
a5556ea5	306
78d7c6d3	307	void AliTrackPoints::MakeITSPoints(AliESDtrack* track)
a5556ea5	308	{
78d7c6d3	309	//Calculates points in ITS
	310	// z=R*Pz/Pt
	311	AliITStrackV2 itstrack(*track,kTRUE);
	312	Double_t x,y,z;
0206ddfb	313	static const Double_t kR[6] = {4.0, 7.0, 14.9, 23.8, 39.1, 43.6};
78d7c6d3	314	for (Int_t i = 0; i < 6; i++)
78d7c6d3	315	{
0206ddfb	316	itstrack.GetGlobalXYZat(kR[i],x,y,z);
78d7c6d3	317	fX[i] = x;
	318	fY[i] = y;
	319	fZ[i] = z;
	320	// Info("MakeITSPoints","X %f Y %f Z %f R asked %f R obtained %f",
	321	// fX[i],fY[i],fZ[i],r[i],TMath::Hypot(fX[i],fY[i]));
	322	}
	323
a5556ea5	324	}
a5556ea5	325
78d7c6d3	326	/***************************************************************/
0206ddfb	327
0f0355ab	328	void AliTrackPoints::MakeITSPointsInnerFromVertexOuterFromTPC(AliESDtrack* track, Float_t mf)
	329	{
	330	//makes trackpoints for ITS
	331	//for 3 inner layers calculates out of the vector at vertex
	332	//for 3 outer ---------------//------------------ at inner TPC
	333
0206ddfb	334	static const Double_t kR[6] = {4.0, 7.0, 14.9, 23.8, 39.1, 43.6};
0f0355ab	335	AliITStrackV2 itstrack(*track,kTRUE);
	336	Double_t x,y,z;
	337	for (Int_t i = 0; i < 3; i++)
	338	{
0206ddfb	339	itstrack.GetGlobalXYZat(kR[i],x,y,z);
0f0355ab	340	fX[i] = x;
	341	fY[i] = y;
	342	fZ[i] = z;
c7ffd78f	343	AliDebug(3,Form("X %f Y %f Z %f R asked %f R obtained %f",
0206ddfb	344	fX[i],fY[i],fZ[i],kR[i],TMath::Hypot(fX[i],fY[i])));
0f0355ab	345	}
	346
	347	for (Int_t i = 3; i < 6; i++)
	348	{
	349	Float_t ax,ay,az;
0206ddfb	350	AliTrackPoints tmptp(1,track,mf,0,kR[i]);
0f0355ab	351	tmptp.PositionAt(0,ax,ay,az);
	352	fX[i] = ax;
	353	fY[i] = ay;
	354	fZ[i] = az;
c7ffd78f	355	AliDebug(3,Form("X %f Y %f Z %f R asked %f R obtained %f",
0206ddfb	356	fX[i],fY[i],fZ[i],kR[i],TMath::Hypot(fX[i],fY[i])));
0f0355ab	357	}
	358
	359	}
	360
0f0355ab	361	/***************************************************************/
0f0355ab	362
a5556ea5	363	void AliTrackPoints::PositionAt(Int_t n, Float_t &x,Float_t &y,Float_t &z)
	364	{
	365	//returns position at point n
5b49024d	366	if ((n<0) \|\| (n>=fN))
a5556ea5	367	{
	368	Error("PositionAt","Point %d out of range",n);
	369	return;
	370	}
	371
	372	x = fX[n];
	373	y = fY[n];
	374	z = fZ[n];
c7ffd78f	375	AliDebug(2,Form("n %d; X %f; Y %f; Z %f",n,x,y,z));
c7ffd78f	376
a5556ea5	377	}
0206ddfb	378
a5556ea5	379	/***************************************************************/
a5556ea5	380
78d7c6d3	381	void AliTrackPoints::Move(Float_t x, Float_t y, Float_t z)
	382	{
	383	//Moves all points about vector
	384	for (Int_t i = 0; i<fN; i++)
	385	{
	386	fX[i]+=x;
	387	fY[i]+=y;
	388	fZ[i]+=z;
	389	}
	390	}
0206ddfb	391
78d7c6d3	392	/***************************************************************/
78d7c6d3	393
a5556ea5	394	Double_t AliTrackPoints::AvarageDistance(const AliTrackPoints& tr)
	395	{
	396	//returns the aritmethic avarage distance between two tracks
	397	// Info("AvarageDistance","Entered");
	398	if ( (fN <= 0) \|\| (tr.fN <=0) )
	399	{
c7ffd78f	400	AliDebug(1,"One of tracks is empty");
a5556ea5	401	return -1;
	402	}
	403
	404	if (fN != tr.fN)
	405	{
	406	Warning("AvarageDistance","Number of points is not equal");
	407	return -1;
	408	}
	409
	410	Double_t sum = 0;
	411	for (Int_t i = 0; i<fN; i++)
	412	{
c7ffd78f	413	AliDebug(10,Form("radii: %f %f",TMath::Hypot(fX[i],fY[i]),TMath::Hypot(tr.fX[i],tr.fY[i])));
a5556ea5	414	// Float_t r1sq = fX[i]fX[i]+fY[i]fY[i];
a5556ea5	415	// Float_t r2sq = tr.fX[i]tr.fX[i]+tr.fY[i]tr.fY[i];
a5556ea5	416
	417
	418	Double_t dx = fX[i]-tr.fX[i];
	419	Double_t dy = fY[i]-tr.fY[i];
	420	Double_t dz = fZ[i]-tr.fZ[i];
	421	sum+=TMath::Sqrt(dxdx + dydy + dz*dz);
	422
c7ffd78f	423	AliDebug(2,Form("Diff: x ,y z: %f , %f, %f",dx,dy,dz));
	424	AliDebug(2,Form("xxyyzz %f %f %f %f %f %f",
	425	fX[i],tr.fX[i],fY[i],tr.fY[i],fZ[i],tr.fZ[i]));
a5556ea5	426	}
	427
	428	Double_t retval = sum/((Double_t)fN);
c7ffd78f	429	AliDebug(1,Form("Avarage distance is %f.",retval));
c7ffd78f	430
a5556ea5	431	return retval;
a5556ea5	432	}
0206ddfb	433
5b49024d	434	/***************************************************************/
	435
	436	void AliTrackPoints::Print(Option_t* /option/) const
	437	{
0206ddfb	438	// Prints the coordinates of the track points
5b49024d	439	Info("Print","There is %d points",fN);
	440	for(Int_t i = 0; i < fN; i++)
	441	{
	442	Info("Print","%d: %f %f %f",i,fX[i],fY[i],fZ[i]);
	443	}
	444
	445	}
	446
a5556ea5	447	/***************************************************************/
a5556ea5	448
78d7c6d3	449	void AliTrackPoints::Testesd(Int_t entr,const char* fname )
a5556ea5	450	{
0206ddfb	451	// This is a former macro converted to member function.
0206ddfb	452	// It is used to test the functionality of the class using ESD tracks
a5556ea5	453	delete gAlice;
	454	gAlice = 0x0;
	455	AliRunLoader* rl = AliRunLoader::Open();
	456	rl->LoadgAlice();
	457
	458	Float_t mf = rl->GetAliRun()->Field()->SolenoidField();
	459
	460
	461	TFile* fFile = TFile::Open(fname);
	462
	463	if (fFile == 0x0)
	464	{
	465	printf("testesd: There is no suche a ESD file\n");
	466	return;
	467	}
	468	AliESD* esd = dynamic_cast<AliESD*>(fFile->Get("0"));
	469	AliESDtrack *t = esd->GetTrack(entr);
	470	if (t == 0x0)
	471	{
	472	::Error("testesd","Can not get track %d",entr);
	473	return;
	474	}
	475
	476
0206ddfb	477	Int_t nTrackPoints = 170;
0206ddfb	478	AliTrackPoints* tp = new AliTrackPoints(nTrackPoints,t,mf,1.);
a5556ea5	479
	480	Float_t xmin = -250;
	481	Float_t xmax = 250;
	482
	483	Float_t ymin = -250;
	484	Float_t ymax = 250;
	485
	486	Float_t zmin = -250;
	487	Float_t zmax = 250;
	488
	489	TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
	490	TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
	491	TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);
	492
	493	TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
	494	TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
	495	TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);
	496
	497	hxyt->SetDirectory(0x0);
	498	hxy->SetDirectory(0x0);
	499	hxyTR->SetDirectory(0x0);
	500
	501	hxzt->SetDirectory(0x0);
	502	hxz->SetDirectory(0x0);
	503	hxzTR->SetDirectory(0x0);
	504
	505	Float_t x,y,z;
	506
0206ddfb	507	for (Int_t i = 0;i<nTrackPoints;i++)
a5556ea5	508	{
	509	Double_t r = 84.1+i;
	510	tp->PositionAt(i,x,y,z);
	511	hxy->Fill(x,y);
	512	hxz->Fill(x,z);
	513	printf("Rdemanded %f\n",r);
	514	printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
	515
	516	}
	517
	518	rl->LoadTrackRefs();
	519	TTree* treeTR = rl->TreeTR();
	520	TBranch* b = treeTR->GetBranch("TPC");
	521
	522	TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
	523	AliTrackReference* tref;
	524	b->SetAddress(&trackrefs);
	525
	526	Int_t tlab = TMath::Abs(t->GetLabel());
	527
	528	Int_t netr = (Int_t)treeTR->GetEntries();
	529	printf("Found %d entries in TR tree\n",netr);
	530
	531	for (Int_t e = 0; e < netr; e++)
	532	{
	533	treeTR->GetEntry(e);
	534	tref = (AliTrackReference*)trackrefs->At(0);
	535	if (tref == 0x0) continue;
	536	if (tref->GetTrack() != tlab) continue;
	537
	538	printf("Found %d entries in TR array\n",trackrefs->GetEntries());
	539
	540	for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
	541	{
	542	tref = (AliTrackReference*)trackrefs->At(i);
	543	if (tref->GetTrack() != tlab) continue;
	544	x = tref->X();
	545	y = tref->Y();
	546	z = tref->Z();
	547	printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
	548
	549	hxzTR->Fill(x,z);
	550	hxyTR->Fill(x,y);
	551	for (Int_t j = 1; j < 10; j++)
	552	{
	553	hxyTR->Fill(x, y+j*0.1);
	554	hxyTR->Fill(x, y-j*0.1);
	555	hxyTR->Fill(x+j*0.1,y);
	556	hxyTR->Fill(x-j*0.1,y);
	557
	558	hxzTR->Fill(x,z-j*0.1);
	559	hxzTR->Fill(x,z+j*0.1);
	560	hxzTR->Fill(x-j*0.1,z);
	561	hxzTR->Fill(x+j*0.1,z);
	562	}
	563	}
	564	break;
	565	}
	566	hxy->Draw("");
	567	// hxzt->Draw("same");
	568	hxyTR->Draw("same");
	569
	570	delete rl;
	571	}
572
a5556ea5	573	/***************************************************************/
a5556ea5	574
78d7c6d3	575	void AliTrackPoints::Testtpc(Int_t entr)
a5556ea5	576	{
0206ddfb	577	// This is a former macro converted to member function.
0206ddfb	578	// It is used to test the functionality of the class using TPC tracks
a5556ea5	579	delete gAlice;
	580	gAlice = 0x0;
	581	AliRunLoader* rl = AliRunLoader::Open();
	582	AliLoader* l = rl->GetLoader("TPCLoader");
	583	rl->LoadgAlice();
7d0f8548	584	AliTracker::SetFieldMap(rl->GetAliRun()->Field(),kTRUE);
a5556ea5	585	l->LoadTracks();
	586	AliTPCtrack* t = new AliTPCtrack();
	587	TBranch* b=l->TreeT()->GetBranch("tracks");
	588	b->SetAddress(&t);
	589	l->TreeT()->GetEntry(entr);
0206ddfb	590	Int_t nTrackPoints = 160;
0206ddfb	591	AliTrackPoints* tp = new AliTrackPoints(nTrackPoints,t,1.);
a5556ea5	592
	593	Float_t xmin = -250;
	594	Float_t xmax = 250;
	595
	596	Float_t ymin = -250;
	597	Float_t ymax = 250;
	598
	599	Float_t zmin = -250;
	600	Float_t zmax = 250;
	601
	602	TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
	603	TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
	604	TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);
	605
	606	TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
	607	TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
	608	TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);
	609
	610	hxyt->SetDirectory(0x0);
	611	hxy->SetDirectory(0x0);
	612	hxyTR->SetDirectory(0x0);
	613
	614	hxzt->SetDirectory(0x0);
	615	hxz->SetDirectory(0x0);
	616	hxzTR->SetDirectory(0x0);
	617
	618	Float_t x,y,z;
	619
0206ddfb	620	for (Int_t i = 0;i<nTrackPoints;i++)
a5556ea5	621	{
	622	Double_t r = 84.1+i;
	623	tp->PositionAt(i,x,y,z);
	624	hxy->Fill(x,y);
	625	hxz->Fill(x,z);
	626	printf("Rdemanded %f\n",r);
	627	printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
5b49024d	628
a5556ea5	629	//BUT they are local!!!!
	630	t->PropagateTo(r);
	631	// Double_t phi = t->Phi();
	632	Double_t rl = TMath::Hypot(t->GetX(),t->GetY());//real radius
	633
	634	Double_t alpha = t->GetAlpha();
	635	Double_t salpha = TMath::Sin(alpha);
	636	Double_t calpha = TMath::Cos(alpha);
	637	x = t->GetX()calpha - t->GetY()salpha;
	638	y = t->GetX()salpha + t->GetY()calpha;
	639	z = t->GetZ();
	640
	641	printf("tx %f ty %f tz %f Rt = %f R from XY %f\n",x,y,z,TMath::Hypot(x,y),rl);
	642
	643	printf("tpz - tz %f\n",z-t->GetZ());
	644	printf("\n");
	645	hxyt->Fill(x,y);
	646	hxzt->Fill(x,z);
	647
	648	}
	649
	650	rl->LoadTrackRefs();
	651	TTree* treeTR = rl->TreeTR();
	652	b = treeTR->GetBranch("TPC");
	653
	654	TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
	655	AliTrackReference* tref;
	656	b->SetAddress(&trackrefs);
	657
	658	Int_t tlab = TMath::Abs(t->GetLabel());
	659
	660	Int_t netr = (Int_t)treeTR->GetEntries();
	661	printf("Found %d entries in TR tree\n",netr);
	662
	663	for (Int_t e = 0; e < netr; e++)
	664	{
	665	treeTR->GetEntry(e);
	666	tref = (AliTrackReference*)trackrefs->At(0);
	667	if (tref == 0x0) continue;
	668	if (tref->GetTrack() != tlab) continue;
	669
	670	printf("Found %d entries in TR array\n",trackrefs->GetEntries());
	671
	672	for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
	673	{
	674	tref = (AliTrackReference*)trackrefs->At(i);
	675	if (tref->GetTrack() != tlab) continue;
	676	x = tref->X();
	677	y = tref->Y();
	678	z = tref->Z();
	679	printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
	680
	681	hxzTR->Fill(x,z);
	682	hxyTR->Fill(x,y);
	683	for (Int_t j = 1; j < 10; j++)
	684	{
	685	hxyTR->Fill(x, y+j*0.1);
	686	hxyTR->Fill(x, y-j*0.1);
	687	hxyTR->Fill(x+j*0.1,y);
	688	hxyTR->Fill(x-j*0.1,y);
	689
	690	hxzTR->Fill(x,z-j*0.1);
	691	hxzTR->Fill(x,z+j*0.1);
	692	hxzTR->Fill(x-j*0.1,z);
693	hxzTR->Fill(x+j*0.1,z);
694	}
695	}
696	break;
697	}
698	hxz->Draw("");
699	// hxzt->Draw("same");
700	hxzTR->Draw("same");
701
702	delete rl;
703	}
78d7c6d3	704