[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 "AliTrackReference.h"

#include "AliTPCtrack.h"
#include "AliESDtrack.h"

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

ClassImp(AliHBTTrackPoints)

Int_t AliHBTTrackPoints::fgDebug = 0;
AliHBTTrackPoints::AliHBTTrackPoints():
 fN(0),
 fX(0x0),
 fY(0x0),
 fZ(0x0)
{
  //constructor
}
AliHBTTrackPoints::AliHBTTrackPoints(const AliHBTTrackPoints& in):
 TObject(in),
 fN(in.fN),
 fX(new Float_t[fN]),
 fY(new Float_t[fN]),
 fZ(new Float_t[fN])
{
//cpy constructor
  for (Int_t i = 0; i < fN; i++)
   {
     fX[i] = in.fX[i];
     fY[i] = in.fY[i];
     fZ[i] = in.fZ[i];
   }
}
/***************************************************************/

AliHBTTrackPoints& AliHBTTrackPoints::operator=(const AliHBTTrackPoints& in)
{
  //Assigment operator
  if(this == &in) return *this;
  TObject::operator=((const TObject&)in);
  
  if (in.fN <= 0)
   {
    delete [] fX;
    delete [] fY;
    delete [] fZ;
    fN = 0;
    fX = 0x0;
    fY = 0x0;
    fZ = 0x0;
    return *this;
   }
   
  if ( fN != in.fN )
   {
    delete [] fX;
    delete [] fY;
    delete [] fZ;
    fN = in.fN;
    fX = new Float_t[fN];
    fY = new Float_t[fN];
    fZ = new Float_t[fN];
   }

  for (Int_t i = 0; i < fN; i++)
   {
     fX[i] = in.fX[i];
     fY[i] = in.fY[i];
     fZ[i] = in.fZ[i];
   }

  return *this;
}
/***************************************************************/

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;
   }
  
  for (Int_t i = 0; i < fN; i++)
   {
     fX[i] = 0.0;
     fY[i] = 0.0;
     fZ[i] = 0.0;
   }
  
  TMath::Hypot(dr,r0)+mf;//just to shut up compilers warning
//  Double_t x;
//  Double_t par[5];
//  track->GetInnerExternalParameters(x,par);     //get properties of the track
//  if (x == 0)
//   {
//     Error("AliHBTTrackPoints","This ESD track does not contain TPC information");
//     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);
}  

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 factorPhi = TMath::ASin( (c2*rc + cst1/rc)/cst2 );//factor phi od rc
     Double_t phi = phi0global + factorPhi - factorPhi0;
     
     Double_t factorZ = TMath::ASin(c2*TMath::Sqrt((rc*rc-dcasq)/cst2))*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]));
      }
   }
}
/***************************************************************/

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

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
  if (fN != tr.fN)
   {
     Error("AvarageDistance","Number of points is not equal");
     return -1;
   }
  if ( (fN <= 0) || (tr.fN <=0) )
   {
     Error("AvarageDistance","One of tracks is empty");
     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()>0)
      {
       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 )
{
//testing routine
  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)
{
//testing routine

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