[u/mrichter/AliRoot.git] / ITS / AliITStrack.cxx

#include <iostream.h>
#include <TMath.h>
#include <TVector.h>
#include <TMatrix.h>
#include <TObjArray.h>

#include "AliITS.h"
#include "AliRun.h"
#include "AliITStrack.h"


ClassImp(AliITStrack)

AliITStrack::AliITStrack() {
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it 
// default constructor   
 
  fvTrack.ResizeTo(5);
  fmCovariance = new TMatrix(5,5);
  flistCluster = new TObjArray; 
  fNumClustInTrack =0;
  fChi2=-1;
  flabel =0; 
  fVertex.ResizeTo(3); 
  fErrorVertex.ResizeTo(3);
  fLayer = -1; 
  ClusterInTrack = new TMatrix(6,9);
  for(Int_t i=0; i<6; i++) (*ClusterInTrack)(i,6)=(*ClusterInTrack)(i,7)=
                           (*ClusterInTrack)(i,8)=-1.;
  rtrack=0.; 
  //alphaprov=-50.; //provvisorio 
  d2.ResizeTo(6);
  tgl2.ResizeTo(6); 
  dtgl.ResizeTo(6);   
}


AliITStrack::AliITStrack(const AliITStrack &cobj) {
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it    

  fvTrack.ResizeTo(5); 
  fmCovariance = new TMatrix(5,5);
  ClusterInTrack = new TMatrix(6,9);
  for(Int_t i=0; i<6; i++) (*ClusterInTrack)(i,6)=(*ClusterInTrack)(i,7)=
                           (*ClusterInTrack)(i,8)=-1.;
  flistCluster = new TObjArray; 
  fVertex.ResizeTo(3); 
  fErrorVertex.ResizeTo(3);
  fVertex = cobj.fVertex;
  fErrorVertex = cobj.fErrorVertex;
  flabel = cobj.flabel;
  fLayer=cobj.fLayer;
  fTPCtrack = cobj.fTPCtrack;
  fNumClustInTrack = cobj.fNumClustInTrack; 
  fChi2= cobj.fChi2;
  fvTrack = cobj.fvTrack;
  rtrack=cobj.rtrack;
  Dv=cobj.Dv;
  Zv=cobj.Zv;
  sigmaDv=cobj.sigmaDv;
  sigmaZv=cobj.sigmaZv;
  d2.ResizeTo(6);
  tgl2.ResizeTo(6); 
  dtgl.ResizeTo(6); 
  d2=cobj.d2;
  tgl2=cobj.tgl2;
  dtgl=cobj.dtgl;
  //alphaprov=cobj.alphaprov;  //provvisorio    
 
  *fmCovariance = *cobj.fmCovariance;
 
  *ClusterInTrack = *cobj.ClusterInTrack;
 
  for(Int_t i=0; i<cobj.flistCluster->GetSize(); i++) 
    flistCluster->AddLast(cobj.flistCluster->At(i));
 
}

AliITStrack::AliITStrack(AliTPCtrack &obj)
{ 
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it 

  fTPCtrack = &obj;
  fvTrack.ResizeTo(5); 
  fVertex.ResizeTo(3); 
  fErrorVertex.ResizeTo(3);

  d2.ResizeTo(6);
  tgl2.ResizeTo(6); 
  dtgl.ResizeTo(6); 
  AliGenerator *gener = gAlice->Generator();
  Float_t Vxg,Vyg,Vzg;
  gener->GetOrigin(Vxg,Vyg,Vzg);
  
 
  fVertex(0)=(Double_t)Vxg;
  fVertex(1)=(Double_t)Vyg; 
  fVertex(2)=(Double_t)Vzg;    
  
  fLayer = 7;
  fmCovariance = new TMatrix(5,5);
  ClusterInTrack = new TMatrix(6,9);
  for(Int_t i=0; i<6; i++) (*ClusterInTrack)(i,6)=(*ClusterInTrack)(i,7)=
                           (*ClusterInTrack)(i,8)=-1.;
  flistCluster = new TObjArray; 
  fNumClustInTrack = 0; 
  LmTPC(); 

}

AliITStrack::~AliITStrack() { 

  //destructor

  if(fmCovariance) delete fmCovariance; 
  if(flistCluster) delete flistCluster; 
  if(ClusterInTrack) delete ClusterInTrack;

}     


void AliITStrack::LmTPC() {

// Transform the TPC state vector from TPC-local to master and build a new state vector ITS-type
// The covariance matrix is also modified accordingly
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it      

  
  TVector tpctrack(5); 
  Double_t Alpha = fTPCtrack->GetAlpha(); 
 
  //printf("LmTPC: Alpha %f\n",Alpha); 

  tpctrack(0) = fTPCtrack->GetY();
  tpctrack(1) = fTPCtrack->GetZ();
  tpctrack(2) = fTPCtrack->GetC();
  tpctrack(3) = fTPCtrack->GetEta();
  tpctrack(4) = fTPCtrack->GetTgl();
  //cout<< " tpctracks = "<<tpctrack(0)<<" "<<tpctrack(1)<<" "<<tpctrack(2)<<" "<<tpctrack(3)<<" "<<tpctrack(4)<<
  //"\n ";

  Double_t xm, ym, zm;
  Double_t sina = TMath::Sin(Alpha);
  Double_t cosa = TMath::Cos(Alpha);
  Double_t xl= fTPCtrack->GetX();
  xm = xl * cosa - tpctrack(0)*sina;
  ym = xl * sina + tpctrack(0)*cosa;
  zm = tpctrack(1);  
  //cout<<" xl e alpha = "<<xl<<" "<<Alpha<<"\n"; getchar(); 

  Double_t x0m,y0m;

  ///////////////////////////////////// determine yo //////////////////////////////////////////////////
  Double_t Vxl=fVertex(0)*cosa+fVertex(1)*sina;
  Double_t Vyl= -fVertex(0)*sina+fVertex(1)*cosa;
  Double_t Xo,Yo, signy;
  Double_t R = 1./tpctrack(2);
  Xo =  tpctrack(3) / tpctrack(2);
  xoTPC=Xo;
  Double_t Yo1, Yo2, diffsq1, diffsq2;  
  Yo1 = tpctrack(0) +  TMath::Sqrt(R*R - (xl-Xo)*(xl-Xo)); 
  Yo2 = tpctrack(0) -  TMath::Sqrt(R*R - (xl-Xo)*(xl-Xo));   
  diffsq1=TMath::Abs((Yo1-Vyl)*(Yo1-Vyl)+(Xo-Vxl)*(Xo-Vxl)-R*R);
  diffsq2=TMath::Abs((Yo2-Vyl)*(Yo2- Vyl)+(Xo-Vxl)*(Xo-Vxl)-R*R);
  if(diffsq1<diffsq2) {Yo=Yo1; signy=1.;} else {Yo=Yo2; signy=-1.;};
  ////////////////////////////////////////////////////////////////////////////////////////////////////

  x0m = Xo * cosa - Yo * sina;
  y0m = Xo * sina + Yo * cosa;  

  rtrack=TMath::Sqrt(xm*xm+ym*ym);
  Double_t pigre=TMath::Pi();
  Double_t phi;
  if(ym == 0. || xm == 0.) {
    if(ym == 0. && xm == 0.) {cout << "  Error in AliITStrack::LmTPC x=0 and y=0 \n"; getchar();}
    if(ym ==0. && xm>0.) phi=0.;
    if(ym==0. && xm<0.) phi=pigre;
    if(xm==0 && ym>0.) phi=pigre/2.;
    if(xm==0 && ym<0.) phi=1.5*pigre;   
  }
  else {
    if (xm>0. && ym>0.) phi=TMath::ATan(ym/xm);
    if (xm<0. && ym>0.) phi=pigre+TMath::ATan(ym/xm);
    if (xm<0. && ym<0.) phi=pigre+TMath::ATan(ym/xm); 
    if (xm>0. && ym<0.) phi=(2*pigre)+TMath::ATan(ym/xm);     
  };
  if(phi<0. || phi>(2*pigre)) {cout<<"attention error on phi in  AliITStrack:LmTPC \n"; getchar();}  

  fvTrack(0)= phi;
  fvTrack(1)=zm;
  fvTrack(3)=tpctrack(4);
  fvTrack(4)=tpctrack(2);
  /*
  //provvisorio
   alphaprov=pigre-PhiDef(x0m,y0m)-PhiDef(fVertex(0),fVertex(1));
  if(alphaprov<0.) alphaprov=alphaprov+2.*pigre;
  //cout<<" PhiDef(x0m,y0m) ="<<PhiDef(x0m,y0m)<<"\n";
  //cout<<" PhiDef(fVertex(0),fVertex(1)) = "<<PhiDef(fVertex(0),fVertex(1))<<"\n";
  //cout<<"alphaprov = "<<alphaprov<<"\n"; getchar();
  //cout<<" fvTrack(0) prima = "<<fvTrack(0)<<"\n";
  fvTrack(0)-=alphaprov;
  if(fvTrack(0)<0.) fvTrack(0)+= 2.*pigre;
  //cout<<" fvTrack(0) dopo = "<<fvTrack(0)<<"  alphaprov = "<<alphaprov<<"\n"; 
  cout<<" fvertex = "<<fVertex(0)<<" "<<fVertex(1)<<" "<<fVertex(2)<<"\n";   
  fVertex(0)=0.;
  fVertex(1)=0.;
/////////////////////////////////////////////////////////////////////////////////////////////////// 
*/ 
  
  Double_t dd=TMath::Sqrt((x0m-fVertex(0))*(x0m-fVertex(0))+(y0m-fVertex(1))*(y0m-fVertex(1)));
  Double_t signdd;
  if (R>0) signdd=1.; else signdd=-1.;
  fvTrack(2)=signdd*dd-R;
  //cout<<" fvertex = "<<fVertex(0)<<" "<<fVertex(1)<<" "<<fVertex(2)<<"\n";
  //cout<<" fvTrack(2) ="<<fvTrack(2)<<"\n"; getchar();
  
  TMatrix localM(5,5);    
  Double_t cov[15];
  fTPCtrack->GetCovariance(cov);
  localM(0,0)=cov[0];
  localM(1,1)=cov[2];
  localM(2,2)=cov[5];
  localM(3,3)=cov[9];
  localM(4,4)=cov[14];
  localM(1,0)=localM(0,1)=cov[1];
  localM(2,0)=localM(0,2)=cov[3];
  localM(2,1)=localM(1,2)=cov[4];
  localM(3,0)=localM(0,3)=cov[6];
  localM(3,1)=localM(1,3)=cov[7];
  localM(3,2)=localM(2,3)=cov[8];
  localM(4,0)=localM(0,4)=cov[10];
  localM(4,1)=localM(1,4)=cov[11];
  localM(4,2)=localM(2,4)=cov[12];
  localM(4,3)=localM(3,4)=cov[13];

  TMatrix F(5,5);
  
  Double_t  dfidy, dDdy, dDdC, dDdeta;

  dfidy=(xm*cosa+ym*sina)/(rtrack*rtrack);
  dDdy=signdd*((y0m-fVertex(1))*cosa-(x0m-fVertex(0))*sina)/dd;
  Double_t dyodr=signy*(R+(xl-Xo)*tpctrack(3))/TMath::Sqrt(R*R-(xl-Xo)*(xl-Xo));
  Double_t dyomdr=sina*tpctrack(3)+cosa*dyodr;
  Double_t dxomdr=cosa*tpctrack(3)-sina*dyodr;
  Double_t ddddR=((x0m-fVertex(0))*dxomdr+(y0m-fVertex(1))*dyomdr)/dd;
  dDdC=-R*R*(signdd*ddddR-1.);
  Double_t dyoldxol=signy*(xl-Xo)/TMath::Sqrt(R*R-(xl-Xo)*(xl-Xo));
  Double_t dxomdeta=R*(cosa-sina*dyoldxol);
  Double_t dyomdeta=R*(sina+cosa*dyoldxol);
  dDdeta=signdd*((x0m-fVertex(0))*dxomdeta+(y0m-fVertex(1))*dyomdeta)/dd;
  F(0,0)=dfidy;
  F(0,1)=F(0,2)=F(0,3)=F(0,4)=0.;
  F(1,1)=1.;
  F(1,0)=F(1,2)=F(1,3)=F(1,4)=0.;
  F(2,0)=dDdy;    
  F(2,2)=dDdC;    
  F(2,3)=dDdeta;  
  F(2,1)=F(2,4)=0.;
  F(3,0)=F(3,1)=F(3,2)=F(3,3)=0.;
  F(3,4)=1.;
  F(4,0)=F(4,1)=F(4,3)=F(4,4)=0.;
  F(4,2)=1.;
// cout<<" Matrice F \n";  
// F.Print(); getchar();
  TMatrix Ft(TMatrix::kTransposed,F);
  TMatrix temp(localM);
  temp*=Ft;
  TMatrix B(F);
  B*=temp;                            // B=F*localM*Ft 
// cout<<" Matrice C\n";  
// B.Print(); getchar();
  *fmCovariance = B;   
}


AliITStrack &AliITStrack::operator=(AliITStrack obj) {
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  

  delete fmCovariance;
  delete flistCluster;
  delete ClusterInTrack;
  fmCovariance = new TMatrix(5,5);
  ClusterInTrack = new TMatrix(6,9);
  flistCluster = new TObjArray; 
  flabel = obj.flabel;
  fTPCtrack = obj.fTPCtrack;
  fvTrack.ResizeTo(5); 
  fNumClustInTrack = obj.fNumClustInTrack; 
  fChi2= obj.fChi2;
  fVertex=obj.fVertex;
  fErrorVertex=obj.fErrorVertex;
  fvTrack = obj.fvTrack; 
  fLayer=obj.fLayer;
  rtrack=obj.rtrack;
  Dv=obj.Dv;
  Zv=obj.Zv;
  sigmaDv=obj.sigmaDv;
  sigmaZv=obj.sigmaZv;
  d2=obj.d2;
  tgl2=obj.tgl2; 
  dtgl=obj.dtgl; 
  //alphaprov=obj.alphaprov; //proviisorio   
    
  *fmCovariance = *obj.fmCovariance;
  *ClusterInTrack = *obj.ClusterInTrack;
  for(Int_t i=0; i<obj.flistCluster->GetSize(); i++) flistCluster->AddLast(obj.flistCluster->At(i));

  return *this;
  
}

void AliITStrack::PutCluster(Int_t layerc, TVector vecclust) {
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
  
  (*ClusterInTrack)(layerc,0) = vecclust(0);
  (*ClusterInTrack)(layerc,1) = vecclust(1);
  (*ClusterInTrack)(layerc,2) = vecclust(2);
  (*ClusterInTrack)(layerc,3) = vecclust(3);
  (*ClusterInTrack)(layerc,4) = vecclust(4);  
  (*ClusterInTrack)(layerc,5) = vecclust(5);
  (*ClusterInTrack)(layerc,6) = vecclust(6);
  (*ClusterInTrack)(layerc,7) = vecclust(7);
  (*ClusterInTrack)(layerc,8) = vecclust(8);

}


void AliITStrack::GetClusters() {

  TMatrix A(*ClusterInTrack);
  TMatrix B(6,3);
  for(Int_t i=0;i<6; i++){
   B(i,0)=A(i,6); B(i,1)=A(i,7); B(i,2)=A(i,8); 
  }
  A.Print();
 // B.Print(); 
 
}


TVector AliITStrack::GetLabTrack(Int_t lay) {
  TVector VecLabel(3);
  VecLabel(0)=( (Float_t) (*ClusterInTrack)(lay,6) );
  VecLabel(1)=( (Float_t) (*ClusterInTrack)(lay,7) );
  VecLabel(2)=( (Float_t) (*ClusterInTrack)(lay,8) );
  return VecLabel;
}

void AliITStrack::Search(TVector VecTotLabref, Long_t &labref, Int_t &freq){
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
// define label

  Int_t vecfreq[18];

  for(Int_t i=0; i<18; i++) vecfreq[i]=0;
  
  for(Int_t i=0; i<18; i++) {
    for(Int_t j=0; j<18; j++) {
      if(VecTotLabref(i) == 0.) VecTotLabref(i)=-3.;  
      if( (VecTotLabref(i)>=0.) && (VecTotLabref(i)==VecTotLabref(j)) ) vecfreq[i]++;    
    }  
  }
  Int_t imax=-1000;
  Long_t  labdefault= (Long_t)1000000.;
  freq=0;
  for(Int_t i=0; i<18; i++) {
    if(vecfreq[i]>freq) {freq=vecfreq[i]; imax=i;}  
  }
  if(imax<0) labref=labdefault; else labref=(Long_t) VecTotLabref(imax);
} 


void AliITStrack::Propagation(Double_t rk) {
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
//Propagation of track 
  Double_t duepi=2.*TMath::Pi();
  Double_t rkm1=rtrack;
//cout<<" rk e rkm1 dentro Propagation   "<<rk<<" "<<rkm1<<"\n";
  
  //
  Double_t Ak=argA(rk), Akm1=argA(rkm1);
  Double_t ak=arga(rk), akm1=arga(rkm1);
  fvTrack(0)+=TMath::ASin(Ak)-TMath::ASin(Akm1);
  if(fvTrack(0)>duepi) fvTrack(0)-=duepi;
  if(fvTrack(0)<0.) fvTrack(0)+=duepi;
	
  Double_t tgl=fvTrack(3);
  Double_t C=fvTrack(4);
  Double_t D=fvTrack(2);
  Double_t Cy=C/2;
  fvTrack(1)+=tgl/Cy*(TMath::ASin(ak)-TMath::ASin(akm1));
  rtrack=rk;
//cout<<"fvTrack(0) fvTrack(1) e rtrack dentro Propagation = "<<fvTrack(0)<<" "<<fvTrack(1)<<" "<<rtrack<<"\n";
//getchar();
	
  Double_t Bk=argB(rk), Bkm1=argB(rkm1);	
  Double_t Ck=argC(rk), Ckm1=argC(rkm1);
  TMatrix F(5,5);
  F(0,2)=Ck/TMath::Sqrt(1.-Ak*Ak) - Ckm1/TMath::Sqrt(1.-Akm1*Akm1);
  F(0,4)=Bk/TMath::Sqrt(1.-Ak*Ak) - Bkm1/TMath::Sqrt(1.-Akm1*Akm1);	
  F(1,2)=tgl*D*(1./rk - 1./rkm1);
  F(1,3) = rk - rkm1;
  F(0,0)=F(1,1)=F(2,2)=F(3,3)=F(4,4)=1.;	
	
  TMatrix Ft(TMatrix::kTransposed,F);
  TMatrix temp(*fmCovariance);
//cout<<" C Matrix prima propagation =\n";
//temp.Print(); getchar();         
  temp*=Ft;
  TMatrix B(F);
  B*=temp;
//cout<<" C Matrix dopo propagation =\n";                            // B=F*C*Ft 
//B.Print(); getchar();
  *fmCovariance = B;    
  	
}
void AliITStrack::AddEL(Double_t signdE, Bool_t flagtot, Double_t mass=0.1396) {
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
//  add energy loss

  TVector s(6);  
  s(0)=0.0026+0.00283; s(1)=0.018; s(2)=0.0094; s(3)=0.0095; s(4)=0.0091; s(5)=0.0087;
//0.00277 is added in the first layer to take into account the energy loss in the beam pipe 
    
  Double_t cl=1.+fvTrack(3)*fvTrack(3);  // cl=1/(cosl)**2 = 1 + (tgl)**2 
  Double_t sqcl=TMath::Sqrt(cl);
  Double_t pt=GetPt();
     
  Double_t p2=pt*pt*cl;
  Double_t E=TMath::Sqrt(p2+mass*mass);
  Double_t beta2=p2/(p2+mass*mass);
  
  Double_t dE;
  if(flagtot) {
    Double_t stot=s(0)+s(1)+s(2)+s(3)+s(4)+s(5);
    dE=0.153/beta2*(log(5940*beta2/(1-beta2)) - beta2)*stot*21.82*sqcl;
  } else {
    dE=0.153/beta2*(log(5940*beta2/(1-beta2)) - beta2)*s(fLayer-1)*21.82*sqcl;
  }   
  dE=signdE*dE/1000.; 
	  
  E+=dE;
  Double_t p=TMath::Sqrt(E*E-mass*mass);   
  Double_t sign=1.;
  if((fvTrack)(4) < 0.) sign=-1.; 
  pt=sign*p/sqcl;  
  Double_t CC=(0.3*0.2)/(pt*100.);
  fvTrack(4)=CC; 
    		
}

void  AliITStrack::Correct(Double_t rk) {
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
// correct track to take into account real geometry detector

  Double_t duepi=2.*TMath::Pi();
  Double_t rkm1=rtrack;
  Double_t Ak=argA(rk), Akm1=argA(rkm1);
  Double_t ak=arga(rk), akm1=arga(rkm1);

  fvTrack(0)+=TMath::ASin(Ak)-TMath::ASin(Akm1);
  if(fvTrack(0)>duepi) fvTrack(0)-=duepi;
  if(fvTrack(0)<0.) fvTrack(0)+=duepi;
	
  Double_t tgl=fvTrack(3);
  Double_t C=fvTrack(4);
  Double_t Cy=C/2;
  fvTrack(1)+=tgl/Cy*(TMath::ASin(ak)-TMath::ASin(akm1));
  rtrack=rk;
  		
}

void AliITStrack::AddMS() {
       
//////////   Modification of the covariance matrix to take into account multiple scattering  ///////////
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it     
   
  TVector s(6);

  s(0)=0.0026+0.00283; s(1)=0.018; s(2)=0.0094; s(3)=0.0095; s(4)=0.0091; s(5)=0.0087;
//0.00277 is added in the first layer to take into account the energy loss in the beam pipe     
  Double_t mass=0.1396;
  Int_t layer=(Int_t)GetLayer();
  
  Double_t tgl=fvTrack(3);
  Double_t cosl=TMath::Cos(TMath::ATan(tgl));	
  Double_t D=fvTrack(2);
  Double_t C=fvTrack(4);
  Double_t Cy=C/2.;
  Double_t Q20=1./(cosl*cosl);	 
  Double_t Q30=C*tgl;
   
  Double_t Q40=Cy*(rtrack*rtrack-D*D)/(1.+ 2.*Cy*D);
  Double_t dd=D+Cy*D*D-Cy*rtrack*rtrack;
  Double_t Q41=-1./cosl*TMath::Sqrt(rtrack*rtrack - dd*dd)/(1.+ 2.*Cy*D);
  
  /*
  Double_t xk=rtrack*TMath::Cos(fvTrack(0));
  Double_t yk=rtrack*TMath::Sin(fvTrack(0));  
  Double_t rvertex=TMath::Sqrt((xk-fVertex(0))*(xk-fVertex(0)) + (yk-fVertex(1))*(yk-fVertex(1)));
  Double_t Q40=Cy*(rvertex*rvertex-D*D)/(1.+ 2.*Cy*D);
  Double_t dd=D+Cy*D*D-Cy*rvertex*rvertex;
  Double_t Q41=-1./cosl*TMath::Sqrt(rvertex*rvertex - dd*dd)/(1.+ 2.*Cy*D); 
  */
   
  TMatrix J(5,2);	
  J(0,0)=0.; J(0,1)=0.;
  J(1,0)=0.; J(1,1)=0.;
  J(2,0)=Q40;    
  J(2,1)=Q41;   	
  J(3,0)=Q20;  J(3,1)=0.;
  J(4,0)=Q30;   J(4,1)=0.;
		 	
  Double_t p2=(GetPt()*GetPt())/(cosl*cosl);
  Double_t beta2=p2/(p2+mass*mass);
  Double_t theta2=14.1*14.1/(beta2*p2*1.e6)*(s(layer-1)/cosl);	
	
  TMatrix Jt(TMatrix::kTransposed,J);	
  TMatrix Q(J,TMatrix::kMult,Jt);
  Q*=theta2;	
		    
  (*fmCovariance)+=Q;
  
}

void AliITStrack::PrimaryTrack() {
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
// calculation of part of covariance matrix for vertex  constraint

  Double_t Rlayer[6];
	
  Rlayer[0]=4.; Rlayer[1]=7.;  Rlayer[2]=14.9;  Rlayer[3]=23.8;  
  Rlayer[4]=39.1;  Rlayer[5]=43.6;
  	
  Double_t Cy=fvTrack(4)/2.;	
  Double_t tgl=(fvTrack(1)-Zv)*Cy/TMath::ASin(Cy*rtrack);
  Double_t rtrack=1.;	
  fvTrack(0)=0.;
  fvTrack(1)=rtrack*tgl; 
  fvTrack(2)=Dv;
  fvTrack(3)=tgl; 
	
  TMatrix  newC(5,5);
  newC(4,4)=(*fmCovariance)(4,4);
  (*fmCovariance)=newC;
  AddEL(1.,1);
  fLayer=0;
  for (Int_t i=0; i<6; i++) {
    Propagation(Rlayer[i]);
    fLayer++;
    d2(i)=(*fmCovariance)(2,2);
    tgl2(i)=(*fmCovariance)(3,3);
    dtgl(i)=(*fmCovariance)(2,3); 
    AddMS();    
    AddEL(-1,0); 	   
  } 		
}	
 
  
Int_t AliITStrack::DoNotCross(Double_t rk) const{
  Double_t C=fvTrack(4);
  Double_t D=fvTrack(2);
  Double_t Cy=C/2.;
  return (TMath::Abs((Cy*rk+(1.+Cy*D)*D/rk)/(1.+2.*Cy*D))>=1.)?1:0;
}
  

Double_t AliITStrack::argA(Double_t rk) const {
  Double_t C=fvTrack(4);
  Double_t D=fvTrack(2);
  Double_t Cy=C/2.;
  Double_t arg=(Cy*rk + (1 + Cy*D)*D/rk)/(1.+ 2.*Cy*D);		
  if (TMath::Abs(arg) < 1.) return arg;
  //cout<<"class AliITSTrack: argA out of range !\n";/* getchar();*/
  return (arg>0) ? 1. : -1.;
}
   
Double_t AliITStrack::arga(Double_t rk) const {
  Double_t C=fvTrack(4);
  Double_t D=fvTrack(2);
  Double_t Cy=C/2.;
  Double_t arg=(rk*rk - D*D)/(1.+ 2.*Cy*D);		
  if (arg<0.) {/*cout<<"class AliITSTrack: arga out of range !\n";*/ arg=0.;} 
  return Cy*TMath::Sqrt(arg);
}
   	
Double_t AliITStrack::argB(Double_t rk) const {
  Double_t C=fvTrack(4);
  Double_t D=fvTrack(2);
  Double_t Cy=C/2.;	   
  return (rk*rk - D*D)/(rk*(1.+ 2.*Cy*D)*(1.+ 2.*Cy*D));
}
   
Double_t AliITStrack::argC(Double_t rk) const {
  Double_t C=fvTrack(4);
  Double_t D=fvTrack(2);
  Double_t Cy=C/2.;		
  return  (1./rk - 2.*Cy*argA(rk)/(1.+ 2.*Cy*D));
}
/*
Double_t AliITStrack::PhiDef(Double_t x, Double_t y){
  Double_t pigre= TMath::Pi();
  Double_t phi;
  if(y == 0. || x == 0.) {
    if(y == 0. && x == 0.) {
      cout << "  Error in AliITStracking::PhiDef x=0 and y=0 \n"; getchar();
    }
    if(y==0. && x>0.) phi=0.;
    if(y==0. && x<0.) phi=pigre;
    if(x==0 && y>0.) phi=pigre/2.;
    if(x==0 && y<0.) phi=1.5*pigre;   
  }
    else {
      if (x>0. && y>0.) phi=TMath::ATan(y/x);
      if (x<0. && y>0.) phi=pigre+TMath::ATan(y/x);
      if (x<0. && y<0.) phi=pigre+TMath::ATan(y/x); 
      if (x>0. && y<0.) phi=(2.*pigre)+TMath::ATan(y/x);     
    }
  if(phi<0. || phi>(2*pigre)) {
    cout<<" Error on phi in  AliITStracking::PhiDef \n"; getchar();
  }  
  return phi;
}
*/
Commit	Line	Data
dce1757b	1	#include <iostream.h>
	2	#include <TMath.h>
	3	#include <TVector.h>
	4	#include <TMatrix.h>
	5	#include <TObjArray.h>
	6
	7	#include "AliITS.h"
	8	#include "AliRun.h"
	9	#include "AliITStrack.h"
	10
	11
	12	ClassImp(AliITStrack)
	13
	14	AliITStrack::AliITStrack() {
	15	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
	16	// default constructor
	17
	18	fvTrack.ResizeTo(5);
	19	fmCovariance = new TMatrix(5,5);
	20	flistCluster = new TObjArray;
	21	fNumClustInTrack =0;
	22	fChi2=-1;
	23	flabel =0;
	24	fVertex.ResizeTo(3);
	25	fErrorVertex.ResizeTo(3);
	26	fLayer = -1;
	27	ClusterInTrack = new TMatrix(6,9);
	28	for(Int_t i=0; i<6; i++) (ClusterInTrack)(i,6)=(ClusterInTrack)(i,7)=
	29	(*ClusterInTrack)(i,8)=-1.;
	30	rtrack=0.;
	31	//alphaprov=-50.; //provvisorio
	32	d2.ResizeTo(6);
	33	tgl2.ResizeTo(6);
	34	dtgl.ResizeTo(6);
	35	}
	36
	37
	38
	39	AliITStrack::AliITStrack(const AliITStrack &cobj) {
	40	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
	41
	42	fvTrack.ResizeTo(5);
	43	fmCovariance = new TMatrix(5,5);
	44	ClusterInTrack = new TMatrix(6,9);
	45	for(Int_t i=0; i<6; i++) (ClusterInTrack)(i,6)=(ClusterInTrack)(i,7)=
	46	(*ClusterInTrack)(i,8)=-1.;
	47	flistCluster = new TObjArray;
	48	fVertex.ResizeTo(3);
	49	fErrorVertex.ResizeTo(3);
	50	fVertex = cobj.fVertex;
	51	fErrorVertex = cobj.fErrorVertex;
	52	flabel = cobj.flabel;
	53	fLayer=cobj.fLayer;
	54	fTPCtrack = cobj.fTPCtrack;
	55	fNumClustInTrack = cobj.fNumClustInTrack;
	56	fChi2= cobj.fChi2;
	57	fvTrack = cobj.fvTrack;
	58	rtrack=cobj.rtrack;
	59	Dv=cobj.Dv;
	60	Zv=cobj.Zv;
	61	sigmaDv=cobj.sigmaDv;
	62	sigmaZv=cobj.sigmaZv;
	63	d2.ResizeTo(6);
	64	tgl2.ResizeTo(6);
65	dtgl.ResizeTo(6);
66	d2=cobj.d2;
67	tgl2=cobj.tgl2;
68	dtgl=cobj.dtgl;
69	//alphaprov=cobj.alphaprov; //provvisorio
70
71	fmCovariance = cobj.fmCovariance;
72
73	ClusterInTrack = cobj.ClusterInTrack;
74
75	for(Int_t i=0; i<cobj.flistCluster->GetSize(); i++)
76	flistCluster->AddLast(cobj.flistCluster->At(i));
77
78	}
79
80	AliITStrack::AliITStrack(AliTPCtrack &obj)
81	{
82	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
83
84	fTPCtrack = &obj;
85	fvTrack.ResizeTo(5);
86	fVertex.ResizeTo(3);
87	fErrorVertex.ResizeTo(3);
88
89	d2.ResizeTo(6);
90	tgl2.ResizeTo(6);
91	dtgl.ResizeTo(6);
92	AliGenerator *gener = gAlice->Generator();
93	Float_t Vxg,Vyg,Vzg;
94	gener->GetOrigin(Vxg,Vyg,Vzg);
95
96
97	fVertex(0)=(Double_t)Vxg;
98	fVertex(1)=(Double_t)Vyg;
99	fVertex(2)=(Double_t)Vzg;
100
101	fLayer = 7;
102	fmCovariance = new TMatrix(5,5);
103	ClusterInTrack = new TMatrix(6,9);
104	for(Int_t i=0; i<6; i++) (ClusterInTrack)(i,6)=(ClusterInTrack)(i,7)=
105	(*ClusterInTrack)(i,8)=-1.;
106	flistCluster = new TObjArray;
107	fNumClustInTrack = 0;
108	LmTPC();
109
110	}
111
112	AliITStrack::~AliITStrack() {
113
114	//destructor
115
116	if(fmCovariance) delete fmCovariance;
117	if(flistCluster) delete flistCluster;
118	if(ClusterInTrack) delete ClusterInTrack;
119
120	}
121
122
123
124	void AliITStrack::LmTPC() {
125
126	// Transform the TPC state vector from TPC-local to master and build a new state vector ITS-type
127	// The covariance matrix is also modified accordingly
128	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
129
130
131	TVector tpctrack(5);
132	Double_t Alpha = fTPCtrack->GetAlpha();
133
134	//printf("LmTPC: Alpha %f\n",Alpha);
135
136	tpctrack(0) = fTPCtrack->GetY();
137	tpctrack(1) = fTPCtrack->GetZ();
138	tpctrack(2) = fTPCtrack->GetC();
139	tpctrack(3) = fTPCtrack->GetEta();
140	tpctrack(4) = fTPCtrack->GetTgl();
141	//cout<< " tpctracks = "<<tpctrack(0)<<" "<<tpctrack(1)<<" "<<tpctrack(2)<<" "<<tpctrack(3)<<" "<<tpctrack(4)<<
142	//"\n ";
143
144	Double_t xm, ym, zm;
145	Double_t sina = TMath::Sin(Alpha);
146	Double_t cosa = TMath::Cos(Alpha);
147	Double_t xl= fTPCtrack->GetX();
148	xm = xl * cosa - tpctrack(0)*sina;
149	ym = xl * sina + tpctrack(0)*cosa;
150	zm = tpctrack(1);
151	//cout<<" xl e alpha = "<<xl<<" "<<Alpha<<"\n"; getchar();
152
153	Double_t x0m,y0m;
154
155	///////////////////////////////////// determine yo //////////////////////////////////////////////////
156	Double_t Vxl=fVertex(0)cosa+fVertex(1)sina;
157	Double_t Vyl= -fVertex(0)sina+fVertex(1)cosa;
158	Double_t Xo,Yo, signy;
159	Double_t R = 1./tpctrack(2);
160	Xo = tpctrack(3) / tpctrack(2);
161	xoTPC=Xo;
162	Double_t Yo1, Yo2, diffsq1, diffsq2;
163	Yo1 = tpctrack(0) + TMath::Sqrt(RR - (xl-Xo)(xl-Xo));
164	Yo2 = tpctrack(0) - TMath::Sqrt(RR - (xl-Xo)(xl-Xo));
165	diffsq1=TMath::Abs((Yo1-Vyl)(Yo1-Vyl)+(Xo-Vxl)(Xo-Vxl)-R*R);
166	diffsq2=TMath::Abs((Yo2-Vyl)(Yo2- Vyl)+(Xo-Vxl)(Xo-Vxl)-R*R);
167	if(diffsq1<diffsq2) {Yo=Yo1; signy=1.;} else {Yo=Yo2; signy=-1.;};
168	////////////////////////////////////////////////////////////////////////////////////////////////////
169
170	x0m = Xo * cosa - Yo * sina;
171	y0m = Xo * sina + Yo * cosa;
172
173	rtrack=TMath::Sqrt(xmxm+ymym);
174	Double_t pigre=TMath::Pi();
175	Double_t phi;
176	if(ym == 0. \|\| xm == 0.) {
177	if(ym == 0. && xm == 0.) {cout << " Error in AliITStrack::LmTPC x=0 and y=0 \n"; getchar();}
178	if(ym ==0. && xm>0.) phi=0.;
179	if(ym==0. && xm<0.) phi=pigre;
180	if(xm==0 && ym>0.) phi=pigre/2.;
181	if(xm==0 && ym<0.) phi=1.5*pigre;
182	}
183	else {
184	if (xm>0. && ym>0.) phi=TMath::ATan(ym/xm);
185	if (xm<0. && ym>0.) phi=pigre+TMath::ATan(ym/xm);
186	if (xm<0. && ym<0.) phi=pigre+TMath::ATan(ym/xm);
187	if (xm>0. && ym<0.) phi=(2*pigre)+TMath::ATan(ym/xm);
188	};
189	if(phi<0. \|\| phi>(2*pigre)) {cout<<"attention error on phi in AliITStrack:LmTPC \n"; getchar();}
190
191	fvTrack(0)= phi;
192	fvTrack(1)=zm;
193	fvTrack(3)=tpctrack(4);
194	fvTrack(4)=tpctrack(2);
195	/*
196	//provvisorio
197	alphaprov=pigre-PhiDef(x0m,y0m)-PhiDef(fVertex(0),fVertex(1));
198	if(alphaprov<0.) alphaprov=alphaprov+2.*pigre;
199	//cout<<" PhiDef(x0m,y0m) ="<<PhiDef(x0m,y0m)<<"\n";
200	//cout<<" PhiDef(fVertex(0),fVertex(1)) = "<<PhiDef(fVertex(0),fVertex(1))<<"\n";
201	//cout<<"alphaprov = "<<alphaprov<<"\n"; getchar();
202	//cout<<" fvTrack(0) prima = "<<fvTrack(0)<<"\n";
203	fvTrack(0)-=alphaprov;
204	if(fvTrack(0)<0.) fvTrack(0)+= 2.*pigre;
205	//cout<<" fvTrack(0) dopo = "<<fvTrack(0)<<" alphaprov = "<<alphaprov<<"\n";
206	cout<<" fvertex = "<<fVertex(0)<<" "<<fVertex(1)<<" "<<fVertex(2)<<"\n";
207	fVertex(0)=0.;
208	fVertex(1)=0.;
209	///////////////////////////////////////////////////////////////////////////////////////////////////
210	*/
211
212	Double_t dd=TMath::Sqrt((x0m-fVertex(0))(x0m-fVertex(0))+(y0m-fVertex(1))(y0m-fVertex(1)));
213	Double_t signdd;
214	if (R>0) signdd=1.; else signdd=-1.;
215	fvTrack(2)=signdd*dd-R;
216	//cout<<" fvertex = "<<fVertex(0)<<" "<<fVertex(1)<<" "<<fVertex(2)<<"\n";
217	//cout<<" fvTrack(2) ="<<fvTrack(2)<<"\n"; getchar();
218
219	TMatrix localM(5,5);
220	Double_t cov[15];
221	fTPCtrack->GetCovariance(cov);
222	localM(0,0)=cov[0];
223	localM(1,1)=cov[2];
224	localM(2,2)=cov[5];
225	localM(3,3)=cov[9];
226	localM(4,4)=cov[14];
227	localM(1,0)=localM(0,1)=cov[1];
228	localM(2,0)=localM(0,2)=cov[3];
229	localM(2,1)=localM(1,2)=cov[4];
230	localM(3,0)=localM(0,3)=cov[6];
231	localM(3,1)=localM(1,3)=cov[7];
232	localM(3,2)=localM(2,3)=cov[8];
233	localM(4,0)=localM(0,4)=cov[10];
234	localM(4,1)=localM(1,4)=cov[11];
235	localM(4,2)=localM(2,4)=cov[12];
236	localM(4,3)=localM(3,4)=cov[13];
237
238	TMatrix F(5,5);
239
240	Double_t dfidy, dDdy, dDdC, dDdeta;
241
242	dfidy=(xmcosa+ymsina)/(rtrack*rtrack);
243	dDdy=signdd((y0m-fVertex(1))cosa-(x0m-fVertex(0))*sina)/dd;
244	Double_t dyodr=signy(R+(xl-Xo)tpctrack(3))/TMath::Sqrt(RR-(xl-Xo)(xl-Xo));
245	Double_t dyomdr=sinatpctrack(3)+cosadyodr;
246	Double_t dxomdr=cosatpctrack(3)-sinadyodr;
247	Double_t ddddR=((x0m-fVertex(0))dxomdr+(y0m-fVertex(1))dyomdr)/dd;
248	dDdC=-RR(signdd*ddddR-1.);
249	Double_t dyoldxol=signy(xl-Xo)/TMath::Sqrt(RR-(xl-Xo)*(xl-Xo));
250	Double_t dxomdeta=R(cosa-sinadyoldxol);
251	Double_t dyomdeta=R(sina+cosadyoldxol);
252	dDdeta=signdd((x0m-fVertex(0))dxomdeta+(y0m-fVertex(1))*dyomdeta)/dd;
253	F(0,0)=dfidy;
254	F(0,1)=F(0,2)=F(0,3)=F(0,4)=0.;
255	F(1,1)=1.;
256	F(1,0)=F(1,2)=F(1,3)=F(1,4)=0.;
257	F(2,0)=dDdy;
258	F(2,2)=dDdC;
259	F(2,3)=dDdeta;
260	F(2,1)=F(2,4)=0.;
261	F(3,0)=F(3,1)=F(3,2)=F(3,3)=0.;
262	F(3,4)=1.;
263	F(4,0)=F(4,1)=F(4,3)=F(4,4)=0.;
264	F(4,2)=1.;
265	// cout<<" Matrice F \n";
266	// F.Print(); getchar();
267	TMatrix Ft(TMatrix::kTransposed,F);
268	TMatrix temp(localM);
269	temp*=Ft;
270	TMatrix B(F);
271	B=temp; // B=FlocalM*Ft
272	// cout<<" Matrice C\n";
273	// B.Print(); getchar();
274	*fmCovariance = B;
275	}
276
277
278	AliITStrack &AliITStrack::operator=(AliITStrack obj) {
279	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
280
281	delete fmCovariance;
282	delete flistCluster;
283	delete ClusterInTrack;
284	fmCovariance = new TMatrix(5,5);
285	ClusterInTrack = new TMatrix(6,9);
286	flistCluster = new TObjArray;
287	flabel = obj.flabel;
288	fTPCtrack = obj.fTPCtrack;
289	fvTrack.ResizeTo(5);
290	fNumClustInTrack = obj.fNumClustInTrack;
291	fChi2= obj.fChi2;
292	fVertex=obj.fVertex;
293	fErrorVertex=obj.fErrorVertex;
294	fvTrack = obj.fvTrack;
295	fLayer=obj.fLayer;
296	rtrack=obj.rtrack;
297	Dv=obj.Dv;
298	Zv=obj.Zv;
299	sigmaDv=obj.sigmaDv;
300	sigmaZv=obj.sigmaZv;
301	d2=obj.d2;
302	tgl2=obj.tgl2;
303	dtgl=obj.dtgl;
304	//alphaprov=obj.alphaprov; //proviisorio
305
306	fmCovariance = obj.fmCovariance;
307	ClusterInTrack = obj.ClusterInTrack;
308	for(Int_t i=0; i<obj.flistCluster->GetSize(); i++) flistCluster->AddLast(obj.flistCluster->At(i));
309
310	return *this;
311
312	}
313
314	void AliITStrack::PutCluster(Int_t layerc, TVector vecclust) {
315	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
316
317	(*ClusterInTrack)(layerc,0) = vecclust(0);
318	(*ClusterInTrack)(layerc,1) = vecclust(1);
319	(*ClusterInTrack)(layerc,2) = vecclust(2);
320	(*ClusterInTrack)(layerc,3) = vecclust(3);
321	(*ClusterInTrack)(layerc,4) = vecclust(4);
322	(*ClusterInTrack)(layerc,5) = vecclust(5);
323	(*ClusterInTrack)(layerc,6) = vecclust(6);
324	(*ClusterInTrack)(layerc,7) = vecclust(7);
325	(*ClusterInTrack)(layerc,8) = vecclust(8);
326
327	}
328
329
330	void AliITStrack::GetClusters() {
331
332	TMatrix A(*ClusterInTrack);
333	TMatrix B(6,3);
334	for(Int_t i=0;i<6; i++){
335	B(i,0)=A(i,6); B(i,1)=A(i,7); B(i,2)=A(i,8);
336	}
337	A.Print();
338	// B.Print();
339
340	}
341
342
343	TVector AliITStrack::GetLabTrack(Int_t lay) {
344	TVector VecLabel(3);
345	VecLabel(0)=( (Float_t) (*ClusterInTrack)(lay,6) );
346	VecLabel(1)=( (Float_t) (*ClusterInTrack)(lay,7) );
347	VecLabel(2)=( (Float_t) (*ClusterInTrack)(lay,8) );
348	return VecLabel;
349	}
350
351	void AliITStrack::Search(TVector VecTotLabref, Long_t &labref, Int_t &freq){
352	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
353	// define label
354
355	Int_t vecfreq[18];
356
357	for(Int_t i=0; i<18; i++) vecfreq[i]=0;
358
359	for(Int_t i=0; i<18; i++) {
360	for(Int_t j=0; j<18; j++) {
361	if(VecTotLabref(i) == 0.) VecTotLabref(i)=-3.;
362	if( (VecTotLabref(i)>=0.) && (VecTotLabref(i)==VecTotLabref(j)) ) vecfreq[i]++;
363	}
364	}
365	Int_t imax=-1000;
366	Long_t labdefault= (Long_t)1000000.;
367	freq=0;
368	for(Int_t i=0; i<18; i++) {
369	if(vecfreq[i]>freq) {freq=vecfreq[i]; imax=i;}
370	}
371	if(imax<0) labref=labdefault; else labref=(Long_t) VecTotLabref(imax);
372	}
373
374
375	void AliITStrack::Propagation(Double_t rk) {
376	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
377	//Propagation of track
378	Double_t duepi=2.*TMath::Pi();
379	Double_t rkm1=rtrack;
380	//cout<<" rk e rkm1 dentro Propagation "<<rk<<" "<<rkm1<<"\n";
381
382	//
383	Double_t Ak=argA(rk), Akm1=argA(rkm1);
384	Double_t ak=arga(rk), akm1=arga(rkm1);
385	fvTrack(0)+=TMath::ASin(Ak)-TMath::ASin(Akm1);
386	if(fvTrack(0)>duepi) fvTrack(0)-=duepi;
387	if(fvTrack(0)<0.) fvTrack(0)+=duepi;
388
389	Double_t tgl=fvTrack(3);
390	Double_t C=fvTrack(4);
391	Double_t D=fvTrack(2);
392	Double_t Cy=C/2;
393	fvTrack(1)+=tgl/Cy*(TMath::ASin(ak)-TMath::ASin(akm1));
394	rtrack=rk;
395	//cout<<"fvTrack(0) fvTrack(1) e rtrack dentro Propagation = "<<fvTrack(0)<<" "<<fvTrack(1)<<" "<<rtrack<<"\n";
396	//getchar();
397
398	Double_t Bk=argB(rk), Bkm1=argB(rkm1);
399	Double_t Ck=argC(rk), Ckm1=argC(rkm1);
400	TMatrix F(5,5);
401	F(0,2)=Ck/TMath::Sqrt(1.-AkAk) - Ckm1/TMath::Sqrt(1.-Akm1Akm1);
402	F(0,4)=Bk/TMath::Sqrt(1.-AkAk) - Bkm1/TMath::Sqrt(1.-Akm1Akm1);
403	F(1,2)=tglD(1./rk - 1./rkm1);
404	F(1,3) = rk - rkm1;
405	F(0,0)=F(1,1)=F(2,2)=F(3,3)=F(4,4)=1.;
406
407	TMatrix Ft(TMatrix::kTransposed,F);
408	TMatrix temp(*fmCovariance);
409	//cout<<" C Matrix prima propagation =\n";
410	//temp.Print(); getchar();
411	temp*=Ft;
412	TMatrix B(F);
413	B*=temp;
414	//cout<<" C Matrix dopo propagation =\n"; // B=FCFt
415	//B.Print(); getchar();
416	*fmCovariance = B;
417
418	}
419	void AliITStrack::AddEL(Double_t signdE, Bool_t flagtot, Double_t mass=0.1396) {
420	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
421	// add energy loss
422
423	TVector s(6);
424	s(0)=0.0026+0.00283; s(1)=0.018; s(2)=0.0094; s(3)=0.0095; s(4)=0.0091; s(5)=0.0087;
425	//0.00277 is added in the first layer to take into account the energy loss in the beam pipe
426
427	Double_t cl=1.+fvTrack(3)fvTrack(3); // cl=1/(cosl)2 = 1 + (tgl)*2
428	Double_t sqcl=TMath::Sqrt(cl);
429	Double_t pt=GetPt();
430
431	Double_t p2=ptptcl;
432	Double_t E=TMath::Sqrt(p2+mass*mass);
433	Double_t beta2=p2/(p2+mass*mass);
434
435	Double_t dE;
436	if(flagtot) {
437	Double_t stot=s(0)+s(1)+s(2)+s(3)+s(4)+s(5);
438	dE=0.153/beta2(log(5940beta2/(1-beta2)) - beta2)stot21.82*sqcl;
439	} else {
440	dE=0.153/beta2(log(5940beta2/(1-beta2)) - beta2)s(fLayer-1)21.82*sqcl;
441	}
442	dE=signdE*dE/1000.;
443
444	E+=dE;
445	Double_t p=TMath::Sqrt(EE-massmass);
446	Double_t sign=1.;
447	if((fvTrack)(4) < 0.) sign=-1.;
448	pt=sign*p/sqcl;
449	Double_t CC=(0.30.2)/(pt100.);
450	fvTrack(4)=CC;
451
452	}
453
454	void AliITStrack::Correct(Double_t rk) {
455	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
456	// correct track to take into account real geometry detector
457
458	Double_t duepi=2.*TMath::Pi();
459	Double_t rkm1=rtrack;
460	Double_t Ak=argA(rk), Akm1=argA(rkm1);
461	Double_t ak=arga(rk), akm1=arga(rkm1);
462
463	fvTrack(0)+=TMath::ASin(Ak)-TMath::ASin(Akm1);
464	if(fvTrack(0)>duepi) fvTrack(0)-=duepi;
465	if(fvTrack(0)<0.) fvTrack(0)+=duepi;
466
467	Double_t tgl=fvTrack(3);
468	Double_t C=fvTrack(4);
469	Double_t Cy=C/2;
470	fvTrack(1)+=tgl/Cy*(TMath::ASin(ak)-TMath::ASin(akm1));
471	rtrack=rk;
472
473	}
474
475	void AliITStrack::AddMS() {
476
477	////////// Modification of the covariance matrix to take into account multiple scattering ///////////
478	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
479
480	TVector s(6);
481
482	s(0)=0.0026+0.00283; s(1)=0.018; s(2)=0.0094; s(3)=0.0095; s(4)=0.0091; s(5)=0.0087;
483	//0.00277 is added in the first layer to take into account the energy loss in the beam pipe
484	Double_t mass=0.1396;
485	Int_t layer=(Int_t)GetLayer();
486
487	Double_t tgl=fvTrack(3);
488	Double_t cosl=TMath::Cos(TMath::ATan(tgl));
489	Double_t D=fvTrack(2);
490	Double_t C=fvTrack(4);
491	Double_t Cy=C/2.;
492	Double_t Q20=1./(cosl*cosl);
493	Double_t Q30=C*tgl;
494
495	Double_t Q40=Cy(rtrackrtrack-DD)/(1.+ 2.Cy*D);
496	Double_t dd=D+CyDD-Cyrtrackrtrack;
497	Double_t Q41=-1./coslTMath::Sqrt(rtrackrtrack - dddd)/(1.+ 2.Cy*D);
498
499	/*
500	Double_t xk=rtrack*TMath::Cos(fvTrack(0));
501	Double_t yk=rtrack*TMath::Sin(fvTrack(0));
502	Double_t rvertex=TMath::Sqrt((xk-fVertex(0))(xk-fVertex(0)) + (yk-fVertex(1))(yk-fVertex(1)));
503	Double_t Q40=Cy(rvertexrvertex-DD)/(1.+ 2.Cy*D);
504	Double_t dd=D+CyDD-Cyrvertexrvertex;
505	Double_t Q41=-1./coslTMath::Sqrt(rvertexrvertex - dddd)/(1.+ 2.Cy*D);
506	*/
507
508	TMatrix J(5,2);
509	J(0,0)=0.; J(0,1)=0.;
510	J(1,0)=0.; J(1,1)=0.;
511	J(2,0)=Q40;
512	J(2,1)=Q41;
513	J(3,0)=Q20; J(3,1)=0.;
514	J(4,0)=Q30; J(4,1)=0.;
515
516	Double_t p2=(GetPt()GetPt())/(coslcosl);
517	Double_t beta2=p2/(p2+mass*mass);
518	Double_t theta2=14.114.1/(beta2p21.e6)(s(layer-1)/cosl);
519
520	TMatrix Jt(TMatrix::kTransposed,J);
521	TMatrix Q(J,TMatrix::kMult,Jt);
522	Q*=theta2;
523
524	(*fmCovariance)+=Q;
525
526	}
527
528	void AliITStrack::PrimaryTrack() {
529	//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
530	// calculation of part of covariance matrix for vertex constraint
531
532	Double_t Rlayer[6];
533
534	Rlayer[0]=4.; Rlayer[1]=7.; Rlayer[2]=14.9; Rlayer[3]=23.8;
535	Rlayer[4]=39.1; Rlayer[5]=43.6;
536
537	Double_t Cy=fvTrack(4)/2.;
538	Double_t tgl=(fvTrack(1)-Zv)Cy/TMath::ASin(Cyrtrack);
539	Double_t rtrack=1.;
540	fvTrack(0)=0.;
541	fvTrack(1)=rtrack*tgl;
542	fvTrack(2)=Dv;
543	fvTrack(3)=tgl;
544
545	TMatrix newC(5,5);
546	newC(4,4)=(*fmCovariance)(4,4);
547	(*fmCovariance)=newC;
548	AddEL(1.,1);
549	fLayer=0;
550	for (Int_t i=0; i<6; i++) {
551	Propagation(Rlayer[i]);
552	fLayer++;
553	d2(i)=(*fmCovariance)(2,2);
554	tgl2(i)=(*fmCovariance)(3,3);
555	dtgl(i)=(*fmCovariance)(2,3);
556	AddMS();
557	AddEL(-1,0);
558	}
559	}
560
561
562	Int_t AliITStrack::DoNotCross(Double_t rk) const{
563	Double_t C=fvTrack(4);
564	Double_t D=fvTrack(2);
565	Double_t Cy=C/2.;
566	return (TMath::Abs((Cyrk+(1.+CyD)D/rk)/(1.+2.Cy*D))>=1.)?1:0;
567	}
568
569
570	Double_t AliITStrack::argA(Double_t rk) const {
571	Double_t C=fvTrack(4);
572	Double_t D=fvTrack(2);
573	Double_t Cy=C/2.;
574	Double_t arg=(Cyrk + (1 + CyD)D/rk)/(1.+ 2.Cy*D);
575	if (TMath::Abs(arg) < 1.) return arg;
576	//cout<<"class AliITSTrack: argA out of range !\n";/* getchar();*/
577	return (arg>0) ? 1. : -1.;
578	}
579
580	Double_t AliITStrack::arga(Double_t rk) const {
581	Double_t C=fvTrack(4);
582	Double_t D=fvTrack(2);
583	Double_t Cy=C/2.;
584	Double_t arg=(rkrk - DD)/(1.+ 2.CyD);
585	if (arg<0.) {/cout<<"class AliITSTrack: arga out of range !\n";/ arg=0.;}
586	return Cy*TMath::Sqrt(arg);
587	}
588
589	Double_t AliITStrack::argB(Double_t rk) const {
590	Double_t C=fvTrack(4);
591	Double_t D=fvTrack(2);
592	Double_t Cy=C/2.;
593	return (rkrk - DD)/(rk(1.+ 2.CyD)(1.+ 2.CyD));
594	}
595
596	Double_t AliITStrack::argC(Double_t rk) const {
597	Double_t C=fvTrack(4);
598	Double_t D=fvTrack(2);
599	Double_t Cy=C/2.;
600	return (1./rk - 2.CyargA(rk)/(1.+ 2.CyD));
601	}
602	/*
603	Double_t AliITStrack::PhiDef(Double_t x, Double_t y){
604	Double_t pigre= TMath::Pi();
605	Double_t phi;
606	if(y == 0. \|\| x == 0.) {
607	if(y == 0. && x == 0.) {
608	cout << " Error in AliITStracking::PhiDef x=0 and y=0 \n"; getchar();
609	}
610	if(y==0. && x>0.) phi=0.;
611	if(y==0. && x<0.) phi=pigre;
612	if(x==0 && y>0.) phi=pigre/2.;
613	if(x==0 && y<0.) phi=1.5*pigre;
614	}
615	else {
616	if (x>0. && y>0.) phi=TMath::ATan(y/x);
617	if (x<0. && y>0.) phi=pigre+TMath::ATan(y/x);
618	if (x<0. && y<0.) phi=pigre+TMath::ATan(y/x);
619	if (x>0. && y<0.) phi=(2.*pigre)+TMath::ATan(y/x);
620	}
621	if(phi<0. \|\| phi>(2*pigre)) {
622	cout<<" Error on phi in AliITStracking::PhiDef \n"; getchar();
623	}
624	return phi;
625	}
626	*/