[u/mrichter/AliRoot.git] / HBTAN / AliHBTTrackPoints.cxx

#include "AliHBTTrackPoints.h"
//_________________________________
////////////////////////////////////////////////////////////
//                                                        //
// class AliHBTTrackPoints                                //
//                                                        //
// 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 <TClonesArray.h>
#include <TFile.h>
#include <TMath.h>

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

ClassImp(AliHBTTrackPoints)

Int_t AliHBTTrackPoints::fgDebug = 0;
AliHBTTrackPoints::AliHBTTrackPoints():
 fN(0),
 fX(0x0),
 fY(0x0),
 fZ(0x0)
{
  //constructor
}
/***************************************************************/

AliHBTTrackPoints::AliHBTTrackPoints(AliHBTTrackPoints::ETypes type, AliESDtrack* track):
 fN(0),
 fX(0x0),
 fY(0x0),
 fZ(0x0)
{
  //constructor
  switch (type)
   {
     case kITS:
       //Up to now we keep only point in pixels
       //Used only in non-id analysis
       fN = 1;
       fX = new Float_t[fN];
       fY = new Float_t[fN];
       fZ = new Float_t[fN];
       MakeITSPoints(track);
       break;

     default:
       Info("AliHBTTrackPoints","Not recognized type");
   }
   
}
/***************************************************************/

AliHBTTrackPoints::AliHBTTrackPoints(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 calculated curvature out of Pt
  //r0 - starting radius
  //dr - calculate points every dr cm, default every 30cm
  if (track == 0x0)
   {
     Error("AliHBTTrackPoints","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
     if (GetDebug() > 3) 
       Warning("AliHBTTrackPoints","This ESD track does not contain TPC information");
       
     fN = 0;
     delete [] fX;
     delete [] fY;
     delete [] fZ;
     fX = fY = fZ = 0x0;
     
     return;
   }
  
  Double_t x;
  Double_t par[5];
  track->GetInnerExternalParameters(x,par);     //get properties of the track
  if (par[4] == 0)
   {
     Error("AliHBTTrackPoints","This ESD track seem not to contain TPC information (curv is 0)");
     return;
   }

  if (mf == 0.0)
   {
     Error("AliHBTTrackPoints","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);
  
}
/***************************************************************/

AliHBTTrackPoints::AliHBTTrackPoints(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("AliHBTTrackPoints","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);
}  
/***************************************************************/

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

void AliHBTTrackPoints::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);
  

  if (GetDebug() > 9) 
    Info("AliHBTTrackPoints","Radius0 %f, Real Radius %f",r0,r); 
  
  if (GetDebug() > 5) 
    Info("AliHBTTrackPoints","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) 
      {
        if (GetDebug() > 1) 
          Warning("AliHBTTrackPoints","ASin argument > 1 %f:",ftmp);
        ftmp=1.0;
      }
     else if (ftmp < -1.0) 
      {
        if (GetDebug() > 1) 
          Warning("AliHBTTrackPoints","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) 
      {
        if (GetDebug() > 1) 
          Warning("AliHBTTrackPoints","Sqrt argument < 0: %f",ftmp);
        ftmp=0.0;
      }
     
     ftmp = c2*TMath::Sqrt(ftmp);
     if (ftmp > 1.0) 
      {
        if (GetDebug() > 1) 
          Warning("AliHBTTrackPoints","ASin argument > 1: %f",ftmp);
        ftmp=1.0;
      }
     else if (ftmp < -1.0) 
      {
        if (GetDebug() > 1) 
          Warning("AliHBTTrackPoints","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);
     
     if ( GetDebug() > 2 )
      {
        Info("AliHBTTrackPoints","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 AliHBTTrackPoints::MakeITSPoints(AliESDtrack* track)
{
//Calculates points in ITS
// z=R*Pz/Pt
 AliITStrackV2 itstrack(*track,kTRUE);
 Double_t x,y,z;
 itstrack.GetGlobalXYZat(4.0,x,y,z);
 fX[0] = x;
 fY[0] = y;
 fZ[0] = z;
 
//  Info("MakeITSPoints","X %f Y %f Z %f R asked %f R obtained %f",
//             fX[0],fY[0],fZ[0],4.0,TMath::Hypot(fX[0],fY[0]));
 
}

/***************************************************************/
void AliHBTTrackPoints::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];
  if ( GetDebug() > 1 )
    {
      Info("AliHBTTrackPoints","n %d; X %f; Y %f; Z %f",n,x,y,z);
    }
}
/***************************************************************/

Double_t AliHBTTrackPoints::AvarageDistance(const AliHBTTrackPoints& tr)
{
  //returns the aritmethic avarage distance between two tracks
//  Info("AvarageDistance","Entered");
  if ( (fN <= 0) || (tr.fN <=0) )
   {
     if (GetDebug()) Warning("AvarageDistance","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++)
   {
     if (GetDebug()>9)
      {
//       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];
       Float_t r1sq = TMath::Hypot(fX[i],fY[i]);
       Float_t r2sq = TMath::Hypot(tr.fX[i],tr.fY[i]);
       Info("AvarageDistance","radii: %f %f",r1sq,r2sq);
      } 

      
     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);
     
     if (GetDebug()>1)
      {
       Info("AvarageDistance","Diff: x ,y z: %f , %f, %f",dx,dy,dz);
       Info("AvarageDistance","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);
  if ( GetDebug() )
    {
      Info("AvarageDistance","Avarage distance is %f.",retval);
    }
  return retval;
}
/***************************************************************/
/***************************************************************/
/***************************************************************/
/***************************************************************/
/***************************************************************/
/***************************************************************/

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


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