[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUTrackerCooked.cxx

//-------------------------------------------------------------------------
//               Implementation of the ITS tracker class
//    The pattern recongintion based on the "cooked covariance" approach
//-------------------------------------------------------------------------

#include <Riostream.h>
#include <TTree.h>
#include <TClonesArray.h>

#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliITSUClusterPix.h"
#include "AliITSUGeomTGeo.h"
#include "AliITSUTrackerCooked.h"
#include "AliITSUTrackCooked.h" 

using std::cout;
using std::endl;

ClassImp(AliITSUTrackerCooked)

//************************************************
// Constants hardcoded for the moment:
//************************************************
// radial positions of layers: default contructor
const 
Double_t klRadius[7]={2.34, 3.15, 3.93, 19.61, 24.55, 34.39, 39.34}; //tdr6
// seed "windows" in z and phi: MakeSeeds
const Double_t kzWin=0.33, kpWin=3.14/4;
// Maximal accepted impact parameters for the seeds 
const Double_t kmaxDCAxy=3.;
const Double_t kmaxDCAz= 3.;

//************************************************
// TODO:
//************************************************
// Seeding:
// Precalculate cylidnrical (r,phi) for the clusters;
// use exact r's for the clusters


AliITSUTrackerCooked::AliITSUlayer
              AliITSUTrackerCooked::fgLayers[AliITSUTrackerCooked::kNLayers];

AliITSUTrackerCooked::AliITSUlayer::AliITSUlayer():
  fR(0),
  fN(0),
  fNsel(0) 
{
  //--------------------------------------------------------------------
  // This default constructor needs to be provided
  //--------------------------------------------------------------------
  for (Int_t i=0; i<kMaxClusterPerLayer; i++) {
    fClusters[i]=0;
    fIndex[i]=-1;
  }
}

void AliITSUTrackerCooked::AliITSUlayer::InsertClusters(TClonesArray *clusters)
{
  //--------------------------------------------------------------------
  // Load clusters to this layer
  //--------------------------------------------------------------------
  Int_t ncl=clusters->GetEntriesFast();
 
  while (ncl--) {
     AliITSUClusterPix *c=(AliITSUClusterPix*)clusters->UncheckedAt(ncl);
     c->GoToFrameGlo();
     //if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
     InsertCluster(new AliITSUClusterPix(*c));
  }

}

void AliITSUTrackerCooked::AliITSUlayer::DeleteClusters()
{
  //--------------------------------------------------------------------
  // Load clusters to this layer
  //--------------------------------------------------------------------
  for (Int_t i=0; i<fN; i++) delete fClusters[i];
  fN=0;
}

Int_t 
AliITSUTrackerCooked::AliITSUlayer::InsertCluster(AliCluster *c) {
  //--------------------------------------------------------------------
  // This function inserts a cluster to this layer in increasing
  // order of the cluster's fZ
  //--------------------------------------------------------------------
  if (fN>=kMaxClusterPerLayer) {
     ::Error("InsertCluster","Too many clusters !\n");
     return 1;
  }
  if (fN==0) fClusters[0]=c;
  else {
     Int_t i=FindClusterIndex(c->GetZ());
     Int_t k=fN-i;
     memmove(fClusters+i+1 ,fClusters+i,k*sizeof(AliCluster*));
     fClusters[i]=c;
  }
  fN++;
  return 0;
}

Int_t 
AliITSUTrackerCooked::AliITSUlayer::FindClusterIndex(Double_t z) const {
  //--------------------------------------------------------------------
  // This function returns the index of the first 
  // with its fZ >= "z". 
  //--------------------------------------------------------------------
  if (fN==0) return 0;

  Int_t b=0;
  if (z <= fClusters[b]->GetZ()) return b;

  Int_t e=b+fN-1;
  if (z > fClusters[e]->GetZ()) return e+1;

  Int_t m=(b+e)/2;
  for (; b<e; m=(b+e)/2) {
    if (z > fClusters[m]->GetZ()) b=m+1;
    else e=m; 
  }
  return m;
}


AliITSUTrackerCooked::AliITSUTrackerCooked(): 
  AliTracker(),fSeeds(0) 
{
  //--------------------------------------------------------------------
  // This default constructor needs to be provided
  //--------------------------------------------------------------------

  AliITSUGeomTGeo *gm  = new AliITSUGeomTGeo(kTRUE,kTRUE);
  AliITSUClusterPix::SetGeom(gm);

  for (Int_t i=0; i<kNLayers; i++) fgLayers[i].SetR(klRadius[i]);

  // Some default primary vertex
  Double_t xyz[]={0.,0.,0.};
  Double_t ers[]={2.,2.,2.};

  SetVertex(xyz,ers);

}

AliITSUTrackerCooked::~AliITSUTrackerCooked() 
{
  //--------------------------------------------------------------------
  // Virtual destructor
  //--------------------------------------------------------------------

  UnloadClusters();

  if (fSeeds) {fSeeds->Delete(); delete fSeeds;} 
  fSeeds=0;
}

static Double_t 
f1(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3)
{
    //-----------------------------------------------------------------
    // Initial approximation of the track curvature
    //-----------------------------------------------------------------
    Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
    Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
                    (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
    Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
                    (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
    
    Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
    
    return -xr*yr/sqrt(xr*xr+yr*yr);
}

static Double_t 
f2(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3)
{
    //-----------------------------------------------------------------
    // Initial approximation of the track curvature times center of curvature
    //-----------------------------------------------------------------
    Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
    Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
                    (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
    Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
                    (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
    
    Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
    
    return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
}

static Double_t 
f3(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t z1, Double_t z2)
{
    //-----------------------------------------------------------------
    // Initial approximation of the tangent of the track dip angle
    //-----------------------------------------------------------------
    return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

Bool_t AliITSUTrackerCooked::
AddCookedSeed(const Float_t r1[3], Int_t l1, Int_t i1, 
              const Float_t r2[3], Int_t l2, Int_t i2,
              const AliCluster *c3,Int_t l3, Int_t i3) 
{
    //--------------------------------------------------------------------
    // This is the main cooking function.
    // Creates seed parameters out of provided clusters.
    //--------------------------------------------------------------------
    Float_t x,a;
    if (!c3->GetXAlphaRefPlane(x,a)) return kFALSE;

    Double_t ca=TMath::Cos(a), sa=TMath::Sin(a);
    Double_t x1 = r1[0]*ca + r1[1]*sa,
             y1 =-r1[0]*sa + r1[1]*ca, z1 = r1[2];
    Double_t x2 = r2[0]*ca + r2[1]*sa,
             y2 =-r2[0]*sa + r2[1]*ca, z2 = r2[2];
    Double_t x3 = x,  y3 = c3->GetY(), z3 = c3->GetZ();

    Double_t par[5];
    par[0]=y3;
    par[1]=z3;
    Double_t crv=f1(x1, y1, x2, y2, x3, y3); //curvature
    Double_t cx0=f2(x1, y1, x2, y2, x3, y3); //curvature*x0
    Double_t tgl12=f3(x1, y1, x2, y2, z1, z2);
    Double_t tgl23=f3(x2, y2, x3, y3, z2, z3);

    Double_t sf=x*crv - cx0;
    if (TMath::Abs(sf) >= kAlmost1) return kFALSE;
    par[2]=sf;

    par[3]=0.5*(tgl12 + tgl23);
    Double_t bz=GetBz();
    par[4]=(TMath::Abs(bz) < kAlmost0Field) ? par[4]=kAlmost0 : crv/(bz*kB2C);

    Double_t cov[15];
    for (Int_t i=0; i<15; i++) cov[i]=0.;
    cov[0] =0.0005*0.0005*10;
    cov[2] =0.0005*0.0005*10;
    cov[5] =0.007*0.007*10;   //FIXME all these lines
    cov[9] =0.007*0.007*10;
    cov[14]=0.1*0.1*10;

    AliITSUTrackCooked *seed=new AliITSUTrackCooked();
    seed->Set(Double_t(x), Double_t(a), par, cov);

    Float_t dz[2]; 
    seed->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
    if (TMath::Abs(dz[0]) > kmaxDCAxy) {delete seed; return kFALSE;} 
    if (TMath::Abs(dz[1]) > kmaxDCAz ) {delete seed; return kFALSE;} 

    seed->SetClusterIndex(l1,i1);
    seed->SetClusterIndex(l2,i2);
    seed->SetClusterIndex(l3,i3);

    fSeeds->AddLast(seed);

    return kTRUE;
}

Int_t AliITSUTrackerCooked::MakeSeeds(Int_t l1, Int_t l2) {
  //--------------------------------------------------------------------
  // This is the main pattern recongition function.
  // Creates seeds out of two clusters and another point.
  //--------------------------------------------------------------------
   Int_t l3=5;

   if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
   fSeeds=new TObjArray(77777);

   Double_t zv=GetZ();

   AliITSUlayer &layer1=fgLayers[l1];
   AliITSUlayer &layer2=fgLayers[l2];
   AliITSUlayer &layer3=fgLayers[l3];
   Double_t r1=layer1.GetR();
   Double_t r2=layer2.GetR();
   Double_t r3=layer3.GetR();

   Int_t nClusters1=layer1.GetNumberOfClusters();
   Int_t nClusters2=layer2.GetNumberOfClusters();
   Int_t nClusters3=layer3.GetNumberOfClusters();
   for (Int_t n1=0; n1<nClusters1; n1++) {
     AliCluster *c1=layer1.GetCluster(n1);
     //
     //Int_t lab=c1->GetLabel(0);
     //
     Double_t z1=c1->GetZ();
     Float_t xyz1[3]; c1->GetGlobalXYZ(xyz1);
     Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
     Double_t zr2=zv + r2/r1*(z1-zv);
     Int_t start2=layer2.FindClusterIndex(zr2-kzWin);
     for (Int_t n2=start2; n2<nClusters2; n2++) {
         AliCluster *c2=layer2.GetCluster(n2);
         //
         //if (c2->GetLabel(0)!=lab) continue;
	 //
         Double_t z2=c2->GetZ();
         
         if (z2 > (zr2+kzWin)) break;  //check in Z
         Float_t xyz2[3]; c2->GetGlobalXYZ(xyz2);
         Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
         if (TMath::Abs(phi2-phi1) > kpWin) continue;  //check in Phi
 
         Double_t zr3=z1 + (r3-r1)/(r2-r1)*(z2-z1);
         Double_t d13=r1-r3, d32=r3-r2, d=r1-r2;
         Double_t phir3=d32/d*phi1 + d13/d*phi2;  // FIXME
         Int_t start3=layer3.FindClusterIndex(zr3-kzWin/2);
         for (Int_t n3=start3; n3<nClusters3; n3++) {
             AliCluster *c3=layer3.GetCluster(n3);
             //
             //if (c3->GetLabel(0)!=lab) continue;
             //
             Double_t z3=c3->GetZ();
         
             if (z3 > (zr3+kzWin/2)) break;  //check in Z
             Float_t xyz3[3]; c3->GetGlobalXYZ(xyz3);
             Double_t phi3=TMath::ATan2(xyz3[1],xyz3[0]);
             if (TMath::Abs(phir3-phi3) > kpWin/100) continue;  //check in Phi

             AliITSUClusterPix cc(*((AliITSUClusterPix*)c2));
             cc.GoToFrameTrk();
             AddCookedSeed(xyz1,l1,n1, xyz3,l3,n3, &cc,l2,n2);

	 }
     }
   }

   return fSeeds->GetEntriesFast();
}

Int_t AliITSUTrackerCooked::Clusters2Tracks(AliESDEvent *event) {
  //--------------------------------------------------------------------
  // This is the main tracking function
  // The clusters must already be loaded
  //--------------------------------------------------------------------
  Int_t l1=6, l2=4;

  // Possibly, create the track "seeds" (combinatorial)
  Int_t nSeeds=MakeSeeds(l1,l2);

  // Possibly, icrement the seeds with additional clusters (Kalman)

  // Possibly, (re)fit the found tracks 


  Int_t nClusters1=fgLayers[l1].GetNumberOfClusters();
  Int_t nClusters2=fgLayers[l2].GetNumberOfClusters();
  cout<<nClusters1<<' '<<nClusters2<<endl;

    Int_t ngood=0;
    for (Int_t s=0; s<nSeeds; s++) {
        AliITSUTrackCooked *track=(AliITSUTrackCooked*)fSeeds->At(s);
        CookLabel(track,0.);
        if (track->GetLabel()>0) {ngood++;
            //cout<<track->GetLabel()<<endl;
        }
        AliESDtrack iotrack;
        iotrack.UpdateTrackParams(track,AliESDtrack::kITSin);
        iotrack.SetLabel(track->GetLabel());
        event->AddTrack(&iotrack);
    }
    cout<<"Good tracks "<<ngood<<endl;
    cout<<"Good tracks/seeds "<<Float_t(ngood)/nSeeds<<endl;

  if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
  fSeeds=0;
    
  return 0;
}

Int_t AliITSUTrackerCooked::PropagateBack(AliESDEvent *event) {
  //--------------------------------------------------------------------
  // Here, we implement the Kalman smoother ?
  // The clusters must already be loaded
  //--------------------------------------------------------------------
    Int_t n=event->GetNumberOfTracks();
    for (Int_t i=0; i<n; i++) {
        AliESDtrack *esdTrack=event->GetTrack(i);
        AliITSUTrackCooked track(*esdTrack);
        esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSout);
    }
    
  return 0;
}

Int_t AliITSUTrackerCooked::RefitInward(AliESDEvent *event) {
  //--------------------------------------------------------------------
  // Some final refit, after the outliers get removed by the smoother ?  
  // The clusters must be loaded
  //--------------------------------------------------------------------
    Int_t n=event->GetNumberOfTracks();
    for (Int_t i=0; i<n; i++) {
        AliESDtrack *esdTrack=event->GetTrack(i);
        AliITSUTrackCooked track(*esdTrack);
        esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSrefit);
    }
    
  return 0;
}

Int_t AliITSUTrackerCooked::LoadClusters(TTree *cTree) {
  //--------------------------------------------------------------------
  // This function reads the ITSU clusters from the tree,
  // sort them, distribute over the internal tracker arrays, etc
  //--------------------------------------------------------------------
  if (!cTree) {
     AliFatal("No cluster tree !");
     return 1;
  }

  //This TClonesArray is not the owner of the clusters
  TClonesArray dummy("AliITSUClusterPix",kMaxClusterPerLayer), *clusters=&dummy;

  for (Int_t i=0; i<kNLayers; i++) {
      TBranch *br = cTree->GetBranch(Form("ITSRecPoints%d",i));
      if (!br) AliFatal(Form("No cluster branch for layer %d",i));
      br->SetAddress(&clusters);
      br->GetEvent(0);
      fgLayers[i].InsertClusters(clusters);
      clusters->Delete();
  }

  return 0;
}

void AliITSUTrackerCooked::UnloadClusters() {
  //--------------------------------------------------------------------
  // This function unloads ITSU clusters from the RAM
  //--------------------------------------------------------------------
  for (Int_t i=0; i<kNLayers; i++) fgLayers[i].DeleteClusters();
}

AliCluster *AliITSUTrackerCooked::GetCluster(Int_t index) const {
  //--------------------------------------------------------------------
  //       Return pointer to a given cluster
  //--------------------------------------------------------------------
    Int_t l=(index & 0xf0000000) >> 28;
    Int_t c=(index & 0x0fffffff) >> 00;
    return fgLayers[l].GetCluster(c);
}
Commit	Line	Data
4fb1e9d1	1	//-------------------------------------------------------------------------
	2	// Implementation of the ITS tracker class
	3	// The pattern recongintion based on the "cooked covariance" approach
	4	//-------------------------------------------------------------------------
	5
	6	#include <Riostream.h>
	7	#include <TTree.h>
	8	#include <TClonesArray.h>
	9
	10	#include "AliLog.h"
	11	#include "AliESDEvent.h"
	12	#include "AliITSUClusterPix.h"
	13	#include "AliITSUGeomTGeo.h"
	14	#include "AliITSUTrackerCooked.h"
	15	#include "AliITSUTrackCooked.h"
	16
d7f21519	17	using std::cout;
	18	using std::endl;
	19
4fb1e9d1	20	ClassImp(AliITSUTrackerCooked)
	21
	22	//************************************************
	23	// Constants hardcoded for the moment:
	24	//************************************************
	25	// radial positions of layers: default contructor
	26	const
	27	Double_t klRadius[7]={2.34, 3.15, 3.93, 19.61, 24.55, 34.39, 39.34}; //tdr6
	28	// seed "windows" in z and phi: MakeSeeds
	29	const Double_t kzWin=0.33, kpWin=3.14/4;
	30	// Maximal accepted impact parameters for the seeds
8cd74e8c	31	const Double_t kmaxDCAxy=3.;
8cd74e8c	32	const Double_t kmaxDCAz= 3.;
4fb1e9d1	33
	34	//************************************************
	35	// TODO:
	36	//************************************************
	37	// Seeding:
	38	// Precalculate cylidnrical (r,phi) for the clusters;
	39	// use exact r's for the clusters
	40
	41
	42	AliITSUTrackerCooked::AliITSUlayer
	43	AliITSUTrackerCooked::fgLayers[AliITSUTrackerCooked::kNLayers];
	44
	45	AliITSUTrackerCooked::AliITSUlayer::AliITSUlayer():
	46	fR(0),
	47	fN(0),
	48	fNsel(0)
	49	{
	50	//--------------------------------------------------------------------
	51	// This default constructor needs to be provided
	52	//--------------------------------------------------------------------
	53	for (Int_t i=0; i<kMaxClusterPerLayer; i++) {
	54	fClusters[i]=0;
	55	fIndex[i]=-1;
	56	}
	57	}
	58
	59	void AliITSUTrackerCooked::AliITSUlayer::InsertClusters(TClonesArray *clusters)
	60	{
	61	//--------------------------------------------------------------------
	62	// Load clusters to this layer
	63	//--------------------------------------------------------------------
	64	Int_t ncl=clusters->GetEntriesFast();
	65
	66	while (ncl--) {
	67	AliITSUClusterPix c=(AliITSUClusterPix)clusters->UncheckedAt(ncl);
	68	c->GoToFrameGlo();
	69	//if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
	70	InsertCluster(new AliITSUClusterPix(*c));
	71	}
	72
	73	}
	74
	75	void AliITSUTrackerCooked::AliITSUlayer::DeleteClusters()
	76	{
	77	//--------------------------------------------------------------------
	78	// Load clusters to this layer
	79	//--------------------------------------------------------------------
	80	for (Int_t i=0; i<fN; i++) delete fClusters[i];
	81	fN=0;
	82	}
	83
	84	Int_t
	85	AliITSUTrackerCooked::AliITSUlayer::InsertCluster(AliCluster *c) {
	86	//--------------------------------------------------------------------
	87	// This function inserts a cluster to this layer in increasing
	88	// order of the cluster's fZ
	89	//--------------------------------------------------------------------
	90	if (fN>=kMaxClusterPerLayer) {
	91	::Error("InsertCluster","Too many clusters !\n");
	92	return 1;
	93	}
	94	if (fN==0) fClusters[0]=c;
	95	else {
	96	Int_t i=FindClusterIndex(c->GetZ());
97	Int_t k=fN-i;
98	memmove(fClusters+i+1 ,fClusters+i,ksizeof(AliCluster));
99	fClusters[i]=c;
100	}
101	fN++;
102	return 0;
103	}
104
105	Int_t
106	AliITSUTrackerCooked::AliITSUlayer::FindClusterIndex(Double_t z) const {
107	//--------------------------------------------------------------------
108	// This function returns the index of the first
109	// with its fZ >= "z".
110	//--------------------------------------------------------------------
111	if (fN==0) return 0;
112
113	Int_t b=0;
114	if (z <= fClusters[b]->GetZ()) return b;
115
116	Int_t e=b+fN-1;
117	if (z > fClusters[e]->GetZ()) return e+1;
118
119	Int_t m=(b+e)/2;
120	for (; b<e; m=(b+e)/2) {
121	if (z > fClusters[m]->GetZ()) b=m+1;
122	else e=m;
123	}
124	return m;
125	}
126
127
128	AliITSUTrackerCooked::AliITSUTrackerCooked():
129	AliTracker(),fSeeds(0)
130	{
131	//--------------------------------------------------------------------
132	// This default constructor needs to be provided
133	//--------------------------------------------------------------------
134
135	AliITSUGeomTGeo *gm = new AliITSUGeomTGeo(kTRUE,kTRUE);
136	AliITSUClusterPix::SetGeom(gm);
137
138	for (Int_t i=0; i<kNLayers; i++) fgLayers[i].SetR(klRadius[i]);
139
140	// Some default primary vertex
141	Double_t xyz[]={0.,0.,0.};
142	Double_t ers[]={2.,2.,2.};
143
144	SetVertex(xyz,ers);
145
146	}
147
148	AliITSUTrackerCooked::~AliITSUTrackerCooked()
149	{
150	//--------------------------------------------------------------------
151	// Virtual destructor
152	//--------------------------------------------------------------------
153
154	UnloadClusters();
155
156	if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
157	fSeeds=0;
158	}
159
160	static Double_t
161	f1(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3)
162	{
163	//-----------------------------------------------------------------
164	// Initial approximation of the track curvature
165	//-----------------------------------------------------------------
166	Double_t d=(x2-x1)(y3-y2)-(x3-x2)(y2-y1);
167	Double_t a=0.5((y3-y2)(y2y2-y1y1+x2x2-x1x1)-
168	(y2-y1)(y3y3-y2y2+x3x3-x2*x2));
169	Double_t b=0.5((x2-x1)(y3y3-y2y2+x3x3-x2x2)-
170	(x3-x2)(y2y2-y1y1+x2x2-x1*x1));
171
172	Double_t xr=TMath::Abs(d/(dx1-a)), yr=d/(dy1-b);
173
174	return -xryr/sqrt(xrxr+yr*yr);
175	}
176
177	static Double_t
178	f2(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3)
179	{
180	//-----------------------------------------------------------------
181	// Initial approximation of the track curvature times center of curvature
182	//-----------------------------------------------------------------
183	Double_t d=(x2-x1)(y3-y2)-(x3-x2)(y2-y1);
184	Double_t a=0.5((y3-y2)(y2y2-y1y1+x2x2-x1x1)-
185	(y2-y1)(y3y3-y2y2+x3x3-x2*x2));
186	Double_t b=0.5((x2-x1)(y3y3-y2y2+x3x3-x2x2)-
187	(x3-x2)(y2y2-y1y1+x2x2-x1*x1));
188
189	Double_t xr=TMath::Abs(d/(dx1-a)), yr=d/(dy1-b);
190
191	return -a/(dy1-b)xr/sqrt(xrxr+yryr);
192	}
193
194	static Double_t
195	f3(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t z1, Double_t z2)
196	{
197	//-----------------------------------------------------------------
198	// Initial approximation of the tangent of the track dip angle
199	//-----------------------------------------------------------------
200	return (z1 - z2)/sqrt((x1-x2)(x1-x2)+(y1-y2)(y1-y2));
201	}
202
203	Bool_t AliITSUTrackerCooked::
204	AddCookedSeed(const Float_t r1[3], Int_t l1, Int_t i1,
205	const Float_t r2[3], Int_t l2, Int_t i2,
206	const AliCluster *c3,Int_t l3, Int_t i3)
207	{
208	//--------------------------------------------------------------------
209	// This is the main cooking function.
210	// Creates seed parameters out of provided clusters.
211	//--------------------------------------------------------------------
212	Float_t x,a;
213	if (!c3->GetXAlphaRefPlane(x,a)) return kFALSE;
214
8cd74e8c	215	Double_t ca=TMath::Cos(a), sa=TMath::Sin(a);
	216	Double_t x1 = r1[0]ca + r1[1]sa,
	217	y1 =-r1[0]sa + r1[1]ca, z1 = r1[2];
	218	Double_t x2 = r2[0]ca + r2[1]sa,
	219	y2 =-r2[0]sa + r2[1]ca, z2 = r2[2];
	220	Double_t x3 = x, y3 = c3->GetY(), z3 = c3->GetZ();
	221
4fb1e9d1	222	Double_t par[5];
8cd74e8c	223	par[0]=y3;
	224	par[1]=z3;
	225	Double_t crv=f1(x1, y1, x2, y2, x3, y3); //curvature
	226	Double_t cx0=f2(x1, y1, x2, y2, x3, y3); //curvature*x0
	227	Double_t tgl12=f3(x1, y1, x2, y2, z1, z2);
	228	Double_t tgl23=f3(x2, y2, x3, y3, z2, z3);
	229
	230	Double_t sf=x*crv - cx0;
	231	if (TMath::Abs(sf) >= kAlmost1) return kFALSE;
	232	par[2]=sf;
	233
4fb1e9d1	234	par[3]=0.5*(tgl12 + tgl23);
	235	Double_t bz=GetBz();
	236	par[4]=(TMath::Abs(bz) < kAlmost0Field) ? par[4]=kAlmost0 : crv/(bz*kB2C);
	237
	238	Double_t cov[15];
	239	for (Int_t i=0; i<15; i++) cov[i]=0.;
	240	cov[0] =0.00050.000510;
	241	cov[2] =0.00050.000510;
	242	cov[5] =0.0070.00710; //FIXME all these lines
	243	cov[9] =0.0070.00710;
	244	cov[14]=0.10.110;
	245
	246	AliITSUTrackCooked *seed=new AliITSUTrackCooked();
	247	seed->Set(Double_t(x), Double_t(a), par, cov);
	248
	249	Float_t dz[2];
	250	seed->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
	251	if (TMath::Abs(dz[0]) > kmaxDCAxy) {delete seed; return kFALSE;}
	252	if (TMath::Abs(dz[1]) > kmaxDCAz ) {delete seed; return kFALSE;}
	253
	254	seed->SetClusterIndex(l1,i1);
	255	seed->SetClusterIndex(l2,i2);
	256	seed->SetClusterIndex(l3,i3);
	257
	258	fSeeds->AddLast(seed);
	259
	260	return kTRUE;
	261	}
	262
	263	Int_t AliITSUTrackerCooked::MakeSeeds(Int_t l1, Int_t l2) {
	264	//--------------------------------------------------------------------
	265	// This is the main pattern recongition function.
	266	// Creates seeds out of two clusters and another point.
	267	//--------------------------------------------------------------------
	268	Int_t l3=5;
	269
	270	if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
	271	fSeeds=new TObjArray(77777);
	272
	273	Double_t zv=GetZ();
	274
	275	AliITSUlayer &layer1=fgLayers[l1];
	276	AliITSUlayer &layer2=fgLayers[l2];
	277	AliITSUlayer &layer3=fgLayers[l3];
	278	Double_t r1=layer1.GetR();
	279	Double_t r2=layer2.GetR();
	280	Double_t r3=layer3.GetR();
	281
	282	Int_t nClusters1=layer1.GetNumberOfClusters();
	283	Int_t nClusters2=layer2.GetNumberOfClusters();
	284	Int_t nClusters3=layer3.GetNumberOfClusters();
	285	for (Int_t n1=0; n1<nClusters1; n1++) {
	286	AliCluster *c1=layer1.GetCluster(n1);
	287	//
	288	//Int_t lab=c1->GetLabel(0);
	289	//
	290	Double_t z1=c1->GetZ();
	291	Float_t xyz1[3]; c1->GetGlobalXYZ(xyz1);
	292	Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
	293	Double_t zr2=zv + r2/r1*(z1-zv);
	294	Int_t start2=layer2.FindClusterIndex(zr2-kzWin);
	295	for (Int_t n2=start2; n2<nClusters2; n2++) {
	296	AliCluster *c2=layer2.GetCluster(n2);
	297	//
298	//if (c2->GetLabel(0)!=lab) continue;
299	//
300	Double_t z2=c2->GetZ();
301
302	if (z2 > (zr2+kzWin)) break; //check in Z
303	Float_t xyz2[3]; c2->GetGlobalXYZ(xyz2);
304	Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
305	if (TMath::Abs(phi2-phi1) > kpWin) continue; //check in Phi
306
307	Double_t zr3=z1 + (r3-r1)/(r2-r1)*(z2-z1);
58c05301	308	Double_t d13=r1-r3, d32=r3-r2, d=r1-r2;
58c05301	309	Double_t phir3=d32/dphi1 + d13/dphi2; // FIXME
4fb1e9d1	310	Int_t start3=layer3.FindClusterIndex(zr3-kzWin/2);
	311	for (Int_t n3=start3; n3<nClusters3; n3++) {
	312	AliCluster *c3=layer3.GetCluster(n3);
	313	//
	314	//if (c3->GetLabel(0)!=lab) continue;
	315	//
	316	Double_t z3=c3->GetZ();
	317
	318	if (z3 > (zr3+kzWin/2)) break; //check in Z
	319	Float_t xyz3[3]; c3->GetGlobalXYZ(xyz3);
	320	Double_t phi3=TMath::ATan2(xyz3[1],xyz3[0]);
58c05301	321	if (TMath::Abs(phir3-phi3) > kpWin/100) continue; //check in Phi
4fb1e9d1	322
	323	AliITSUClusterPix cc(((AliITSUClusterPix)c2));
	324	cc.GoToFrameTrk();
	325	AddCookedSeed(xyz1,l1,n1, xyz3,l3,n3, &cc,l2,n2);
	326
	327	}
	328	}
	329	}
	330
	331	return fSeeds->GetEntriesFast();
	332	}
	333
	334	Int_t AliITSUTrackerCooked::Clusters2Tracks(AliESDEvent *event) {
	335	//--------------------------------------------------------------------
	336	// This is the main tracking function
	337	// The clusters must already be loaded
	338	//--------------------------------------------------------------------
	339	Int_t l1=6, l2=4;
	340
	341	// Possibly, create the track "seeds" (combinatorial)
	342	Int_t nSeeds=MakeSeeds(l1,l2);
	343
	344	// Possibly, icrement the seeds with additional clusters (Kalman)
	345
	346	// Possibly, (re)fit the found tracks
	347
	348
	349	Int_t nClusters1=fgLayers[l1].GetNumberOfClusters();
	350	Int_t nClusters2=fgLayers[l2].GetNumberOfClusters();
	351	cout<<nClusters1<<' '<<nClusters2<<endl;
	352
	353	Int_t ngood=0;
	354	for (Int_t s=0; s<nSeeds; s++) {
	355	AliITSUTrackCooked track=(AliITSUTrackCooked)fSeeds->At(s);
	356	CookLabel(track,0.);
	357	if (track->GetLabel()>0) {ngood++;
	358	//cout<<track->GetLabel()<<endl;
	359	}
	360	AliESDtrack iotrack;
	361	iotrack.UpdateTrackParams(track,AliESDtrack::kITSin);
	362	iotrack.SetLabel(track->GetLabel());
	363	event->AddTrack(&iotrack);
	364	}
	365	cout<<"Good tracks "<<ngood<<endl;
	366	cout<<"Good tracks/seeds "<<Float_t(ngood)/nSeeds<<endl;
	367
	368	if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
	369	fSeeds=0;
	370
	371	return 0;
	372	}
	373
	374	Int_t AliITSUTrackerCooked::PropagateBack(AliESDEvent *event) {
	375	//--------------------------------------------------------------------
	376	// Here, we implement the Kalman smoother ?
	377	// The clusters must already be loaded
	378	//--------------------------------------------------------------------
	379	Int_t n=event->GetNumberOfTracks();
	380	for (Int_t i=0; i<n; i++) {
	381	AliESDtrack *esdTrack=event->GetTrack(i);
	382	AliITSUTrackCooked track(*esdTrack);
	383	esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSout);
	384	}
	385
386	return 0;
387	}
388
389	Int_t AliITSUTrackerCooked::RefitInward(AliESDEvent *event) {
390	//--------------------------------------------------------------------
391	// Some final refit, after the outliers get removed by the smoother ?
392	// The clusters must be loaded
393	//--------------------------------------------------------------------
394	Int_t n=event->GetNumberOfTracks();
395	for (Int_t i=0; i<n; i++) {
396	AliESDtrack *esdTrack=event->GetTrack(i);
397	AliITSUTrackCooked track(*esdTrack);
398	esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSrefit);
399	}
400
401	return 0;
402	}
403
404	Int_t AliITSUTrackerCooked::LoadClusters(TTree *cTree) {
405	//--------------------------------------------------------------------
406	// This function reads the ITSU clusters from the tree,
407	// sort them, distribute over the internal tracker arrays, etc
408	//--------------------------------------------------------------------
409	if (!cTree) {
410	AliFatal("No cluster tree !");
411	return 1;
412	}
413
414	//This TClonesArray is not the owner of the clusters
415	TClonesArray dummy("AliITSUClusterPix",kMaxClusterPerLayer), *clusters=&dummy;
416
417	for (Int_t i=0; i<kNLayers; i++) {
418	TBranch *br = cTree->GetBranch(Form("ITSRecPoints%d",i));
419	if (!br) AliFatal(Form("No cluster branch for layer %d",i));
420	br->SetAddress(&clusters);
421	br->GetEvent(0);
422	fgLayers[i].InsertClusters(clusters);
423	clusters->Delete();
424	}
425
426	return 0;
427	}
428
429	void AliITSUTrackerCooked::UnloadClusters() {
430	//--------------------------------------------------------------------
431	// This function unloads ITSU clusters from the RAM
432	//--------------------------------------------------------------------
433	for (Int_t i=0; i<kNLayers; i++) fgLayers[i].DeleteClusters();
434	}
435
436	AliCluster *AliITSUTrackerCooked::GetCluster(Int_t index) const {
437	//--------------------------------------------------------------------
438	// Return pointer to a given cluster
439	//--------------------------------------------------------------------
440	Int_t l=(index & 0xf0000000) >> 28;
441	Int_t c=(index & 0x0fffffff) >> 00;
442	return fgLayers[l].GetCluster(c);
443	}