--- /dev/null
+////////////////////////////////////////////////
+// Reconstructed point class for set:ITS //
+////////////////////////////////////////////////
+
+#include <TMatrix.h>
+
+#include "AliITSiotrack.h"
+
+ClassImp(AliITSiotrack)
+
+AliITSiotrack::AliITSiotrack() {
+
+ // default creator
+ fLab=-3;
+ fX=fZ=fY=0.;
+ fPx=fPy=fPz=0.;
+ for (Int_t i=0;i<6;i++) {fIdModules[i]=fIdPoints[i]=-1; fIdPoints[i]=-1;}
+ fStateVPhi=0.; fStateVZ=0.; fStateVD=0.; fStateVTgl=0.; fStateVC=0.;
+ fRadius=0.; fCharge=0;
+ fC00=fC10=fC11=fC20=fC21=fC22=fC30=fC31=fC32=fC33=fC40=fC41=fC42=fC43=fC44=0.;
+
+}
+
+
+void AliITSiotrack::SetCovMatrix(TMatrix *cov) {
+
+ fC00=(*cov)(0,0);
+ fC10=(*cov)(1,0);
+ fC11=(*cov)(1,1);
+ fC20=(*cov)(2,0);
+ fC21=(*cov)(2,1);
+ fC22=(*cov)(2,2);
+ fC30=(*cov)(3,0);
+ fC31=(*cov)(3,1);
+ fC32=(*cov)(3,2);
+ fC33=(*cov)(3,3);
+ fC40=(*cov)(4,0);
+ fC41=(*cov)(4,1);
+ fC42=(*cov)(4,2);
+ fC43=(*cov)(4,3);
+ fC44=(*cov)(4,4);
+
+}
+
+
+Double_t * AliITSiotrack::GetCovMatrix() {
+
+ Double_t covar[25];
+
+ covar[0]=fC00;
+ covar[1]=fC10;
+ covar[2]=fC20;
+ covar[3]=fC30;
+ covar[4]=fC40;
+ covar[5]=fC10;
+ covar[6]=fC11;
+ covar[7]=fC21;
+ covar[8]=fC31;
+ covar[9]=fC41;
+ covar[10]=fC20;
+ covar[11]=fC21;
+ covar[12]=fC22;
+ covar[13]=fC32;
+ covar[14]=fC42;
+ covar[15]=fC30;
+ covar[16]=fC31;
+ covar[17]=fC32;
+ covar[18]=fC33;
+ covar[19]=fC43;
+ covar[20]=fC40;
+ covar[21]=fC41;
+ covar[22]=fC42;
+ covar[23]=fC43;
+ covar[24]=fC44;
+
+ return covar;
+
+}
+/*
+Double_t * AliITSiotrack::GetStateVector() {
+
+ Double_t statev[5];
+
+ statev[0]=fStateVPhi;
+ statev[1]=fStateVZ;
+ statev[2]=fStateVD;
+ statev[3]=fStateVTgl;
+ statev[4]=fStateVC;
+
+ return statev;
+
+} */
+
+
--- /dev/null
+#ifndef ALIITSIOTRACK_H
+#define ALIITSIOTRACK_H
+
+////////////////////////////////////////////////////
+// Reconstructed space point class for set:ITS //
+////////////////////////////////////////////////////
+
+#include <TObject.h>
+
+class AliITSiotrack : public TObject {
+
+
+ public:
+
+ AliITSiotrack(); // constructor
+ virtual ~AliITSiotrack() {}; // distructor
+ Int_t GetLabel() const {return fLab;} // get track label
+ Int_t GetIdPoint(Int_t i) const {return fIdPoints[i];} // get the identification number for the point
+ Int_t GetIdModule(Int_t i) const {return fIdModules[i];} // get the module number for the point
+ //Double_t * GetStateVector(); // get the state vector (5 elements)
+ Double_t * GetCovMatrix(); // get the covariance matrix as 15 elements array
+ Double_t GetStatePhi() const {return fStateVPhi;} // gets the Phi angle of the state vector
+ Double_t GetStateZ() const {return fStateVZ;} // gets the Z cohordinate of the state vector
+ Double_t GetStateD() const {return fStateVD;} // gets the radial impact parameter of the state vector
+ Double_t GetStateC() const {return fStateVC;} // gets the curvature C of the state vector
+ Double_t GetStateTgl() const {return fStateVTgl;} // gets the dip angle tangent of the state vector
+ Double_t GetRadius() const {return fRadius;} // gets the radius corresponding to the state vector
+ Int_t GetCharge() const {return fCharge;} // gets the particle charge
+ Float_t GetX() const {return fX;} // gets the x cohordinate of the found vertex
+ Float_t GetZ() const {return fZ;} // gets the z cohordinate of the found vertex
+ Float_t GetY() const {return fY;} // gets the y cohordinate of the found vertex
+ Float_t GetPx() const {return fPx;} // gets the x momentum component at the found vertex
+ Float_t GetPy() const {return fPy;} // gets the y momentum component at the found vertex
+ Float_t GetPz()const {return fPz;} // gets the z momentum component at the found vertex
+
+ void SetCovMatrix(TMatrix *cov); // sets the covariance matrix
+ void SetLabel(Int_t lab) {fLab=lab;} // sets the track label
+ void SetIdPoint(Int_t i,Int_t pnt) {fIdPoints[i]=pnt;} // set the identification number for the point
+ void SetIdModule(Int_t i,Int_t mod) {fIdModules[i]=mod;} // set the module number for the point
+
+ void SetStatePhi(Double_t phi) {fStateVPhi=phi;} // sets the Phi angle of the state vector
+ void SetStateZ(Double_t z) {fStateVZ=z;} // sets the Z cohordinate of the state vector
+ void SetStateD(Double_t d) {fStateVD=d;} // sets the radial impact parameter of the state vector
+ void SetStateTgl(Double_t tgl) {fStateVTgl=tgl;} // sets the dip angle tangent of the state vector
+ void SetStateC(Double_t c) {fStateVC=c;} // sets the curvature C of the state vector
+ void SetRadius(Double_t r) {fRadius= r;} // sets the radius corresponding to the state vector
+ void SetCharge(Int_t charge) {fCharge=charge;} // sets the particle charge
+
+ void SetX(Float_t x){fX=x;} // sets the x cohordinate of the found vertex
+ void SetZ(Float_t z){fZ=z;} // sets the z cohordinate of the found vertex
+ void SetY(Float_t y){fY=y;} // sets the y cohordinate of the found vertex
+ void SetPx(Float_t px) {fPx=px;} // sets the x momentum component at the found vertex
+ void SetPy(Float_t py) {fPy=py;} // sets the y momentum component at the found vertex
+ void SetPz(Float_t pz) {fPz=pz;} // sets the z momentum component at the found vertex
+
+ private:
+
+ Int_t fLab; // label of reconstructed track
+ Float_t fX ; // x cohordinate of the found vertex
+ Float_t fY ; // y cohordinate of the found vertex
+ Float_t fZ ; // z cohordinate of the found vertex
+ Float_t fPx; // x component of track momentum at the found vertex
+ Float_t fPy; // y component of track momentum at the found vertex
+ Float_t fPz; // z component of track momentum at the found vertex
+
+ //
+ Int_t fIdPoints[6]; // points assigned to the track (entry # in fRecPoints is given by module #)
+ Int_t fIdModules[6]; // module # corresponding to each point belonging to the track
+
+ Double_t fStateVPhi; //
+ Double_t fStateVZ; // state vector components
+ Double_t fStateVD; //
+ Double_t fStateVTgl; //
+ Double_t fStateVC; //
+
+ Double_t fRadius; // distance of the point from the origin
+ Int_t fCharge; // particle charge
+
+// Covariance matrix
+ Double_t fC00;
+ Double_t fC10, fC11;
+ Double_t fC20, fC21, fC22;
+ Double_t fC30, fC31, fC32, fC33;
+ Double_t fC40, fC41, fC42, fC43, fC44;
+
+ ClassDef(AliITSiotrack,1) // AliITSiotrack class
+};
+
+#endif
+
+
+
--- /dev/null
+#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;
+}
+*/
--- /dev/null
+#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;
+}
+*/
--- /dev/null
+#ifndef ALIITSTRACK_H
+#define ALIITSTRACK_H
+
+#include <TObject.h>
+#include <TMatrix.h>
+#include <TVector.h>
+
+#include "../TPC/AliTPCtrack.h"
+
+class TObjArray;
+// ITS Track Class
+//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
+//
+class AliITStrack : public TObject {
+
+public:
+
+ AliITStrack() ;
+ AliITStrack(AliTPCtrack &obj);
+ AliITStrack(const AliITStrack &cobj);
+ AliITStrack &operator=(AliITStrack obj);
+ ~AliITStrack();
+ Float_t &operator()(Int_t i) const { return fvTrack(i);}
+ Float_t &operator()(Int_t r, Int_t c) const {return (*fmCovariance)(r,c);}
+ Int_t GetNumClust() { return fNumClustInTrack;}
+ void AddClustInTrack() { fNumClustInTrack++;}
+ TObjArray *GetListOfCluster() { return flistCluster;}
+ void SetChi2(Double_t chi2) { fChi2 = chi2;}
+ Double_t GetChi2() { return fChi2;}
+ Double_t GetZ() const {return fvTrack(1);}
+ Double_t GetTgl() const {return fvTrack(3);}
+ Double_t Getrtrack() const{return rtrack;}
+ Double_t Getphi() const{return fvTrack(0);}
+ Double_t GetC() const {return fvTrack(4);}
+ Double_t GetD() const{return fvTrack(2);}
+ Double_t GetPt() const {return 0.3*0.2/(fvTrack(4)*100.);}
+ void SetVertex(TVector &vert) { for(Int_t i=0;i<3;i++) fVertex(i) = vert(i);}
+ void SetErrorVertex(TVector &evert) {for(Int_t i=0;i<3;i++) fErrorVertex(i) = evert(i);}
+
+ void LmTPC(); // trasform state vector and covariance matrix from local TPC to master
+ TVector GetVertex() { return fVertex;}
+ TVector GetErrorVertex() { return fErrorVertex;}
+ Long_t GetLabel() { return flabel;}
+ void SetLabel(Long_t label) { flabel = label;}
+ Int_t GetLayer() { return fLayer;}
+ TMatrix &GetCMatrix() { return *fmCovariance;}
+ TVector &GetVector() { return fvTrack;}
+ void SetLayer(Int_t layer) { fLayer = layer;}
+ AliTPCtrack *GetTPCtrack() { return fTPCtrack;}
+
+ void PutCluster(Int_t layerc, TVector vecclust);
+ TVector GetLabTrack(Int_t lay);
+ void Search(TVector VecTotLabref, Long_t &labref, Int_t &freq);
+ Float_t GetZclusterTrack(Int_t lay) {return ((Float_t) (*ClusterInTrack)(lay,2));}
+ void GetClusters();
+ Int_t GetLabTPC() {return (*fTPCtrack).GetLabel();}
+ Int_t GetIdPoint(Int_t lay) {return ((Int_t) (*ClusterInTrack)(lay,4));}
+ Int_t GetIdModule(Int_t lay) {return ((Int_t) (*ClusterInTrack)(lay,5));}
+
+
+ Int_t DoNotCross(Double_t rk) const;
+ Double_t argA(Double_t rk) const;
+ Double_t arga(Double_t rk) const;
+ Double_t argB(Double_t rk) const;
+ Double_t argC(Double_t rk) const;
+ void Propagation(Double_t rk) ;
+
+ Double_t GetSigmaphi() const{return ((Double_t)(*fmCovariance)(0,0));}
+ Double_t GetSigmaZ() const{return ((Double_t)(*fmCovariance)(1,1));}
+ void AddEL(Double_t signdE, Bool_t flagtot, Double_t mass=0.1396);
+ void AddMS();
+ void Correct(Double_t rk);
+ void SetDv(Double_t x) {Dv=x;}
+ void SetZv(Double_t x) {Zv=x;}
+ Double_t GetDv() {return Dv;}
+ Double_t GetZv() {return Zv;}
+ void SetsigmaDv( Double_t x) {sigmaDv=x;}
+ void SetsigmaZv( Double_t x) {sigmaZv=x;}
+ Double_t GetsigmaDv() {return sigmaDv;}
+ Double_t GetsigmaZv() {return sigmaZv;}
+ void PrimaryTrack();
+ void Setd2(TVector &x) {for(Int_t i=0; i<6; i++){d2(i)=x(i);}}
+ void Settgl2(TVector &x) {for(Int_t i=0; i<6; i++){tgl2(i)=x(i);}}
+ void Setdtgl(TVector &x) {for(Int_t i=0; i<6; i++){dtgl(i)=x(i);}}
+ TVector Getd2() { return d2;}
+ TVector Gettgl2() { return tgl2;}
+ TVector Getdtgl() { return dtgl;}
+ Double_t Getd2(Int_t i){return (Double_t)d2(i);}
+ Double_t Gettgl2(Int_t i){return (Double_t)tgl2(i);}
+ Double_t Getdtgl(Int_t i){return (Double_t)dtgl(i);}
+ Double_t GetxoTPC() {return xoTPC;}
+ // Double_t PhiDef(Double_t x, Double_t y);
+ // Double_t Getalphaprov(){return alphaprov;} //provvirio
+////////////////////////////////////////////////////////////////////////////////////////
+
+ private:
+
+ AliTPCtrack *fTPCtrack; // reference to TPC track
+
+ TVector fvTrack; // state vector: |phi/z/D/tgl/C
+ Double_t rtrack; // radius of courrent layer
+ TMatrix *fmCovariance; // Covariance Matrix
+
+ Double_t fChi2; // fChi^2 of track
+ TObjArray *flistCluster; // list of clusters of the track
+ Int_t fNumClustInTrack; // total number of clusters
+ Long_t flabel; // label of the track
+ TVector fVertex; // vertex coordinates of the track
+ TVector fErrorVertex; // error on the vertex coordinates
+ Int_t fLayer; // current Layer of the track
+ TMatrix *ClusterInTrack; // matrix of clusters belonging to the track
+ // row index = layer-1;
+ // cols index = master coordinates of the clusters
+
+
+ Double_t Dv; // radial impact parameter for vertex constraint
+ Double_t Zv; // longitudinal impact parameter for vertex constraint
+ Double_t sigmaDv; // sigma for Dv extraction
+ Double_t sigmaZv; // sigma for Zv extraction
+ TVector d2; // C(2,2) per primary track
+ TVector tgl2; // C(3,3) per primary track
+ TVector dtgl; // C(2,3) per primary track
+
+ Double_t xoTPC;
+
+ //Double_t alphaprov; //provviosorio
+
+
+ ClassDef(AliITStrack, 1)
+
+};
+
+#endif
+
--- /dev/null
+#include <TMath.h>
+#include <TList.h>
+#include <TTree.h>
+#include <TVector.h>
+#include <TMatrix.h>
+#include <TObjectTable.h>
+
+#include "AliITStracking.h"
+#include "AliRun.h"
+#include "AliITS.h"
+#include "AliITSgeom.h"
+#include "AliITSRecPoint.h"
+#include "AliITStrack.h"
+
+ClassImp(AliITStracking)
+
+
+AliITStracking::AliITStracking(TList *trackITSlist, AliITStrack *reference,
+ AliITS *aliITS, TObjArray *rpoints, Double_t Ptref, Int_t **vettid, Bool_t flagvert) {
+/////////////////////// This function perform the tracking in ITS detectors /////////////////////
+/////////////////////// reference is a pointer to the final best track /////////////////////
+//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
+// The authors thank Mariana Bondila to have help them to resolve some problems. July-2000
+
+
+ Rlayer[0]=4.; Rlayer[1]=7.; Rlayer[2]=14.9; Rlayer[3]=23.8; Rlayer[4]=39.1; Rlayer[5]=43.6;
+
+ for(Int_t index =0; index<trackITSlist->GetSize(); index++) {
+ AliITStrack *trackITS = (AliITStrack *) trackITSlist->At(index);
+ if((*trackITS).GetLayer()==7) reference->SetChi2(10.223e140);
+ //cout <<" Layer inizio = "<<(*trackITS).GetLayer()<<"\n";
+ // cout<<"fvtrack =" <<"\n";
+ // cout << (*trackITS)(0) << " "<<(*trackITS)(1)<<" "<<(*trackITS)(2)<<" "<<(*trackITS)(3)<<" "<<(*trackITS)(4)<<"\n";
+ // cout<< " rtrack = "<<(*trackITS).Getrtrack()<<"\n";
+ // getchar();
+ Double_t Chi2Now, Chi2Ref;
+ if((*trackITS).GetLayer()==1 ) {
+ Chi2Now = trackITS->GetChi2();
+ Float_t NumClustNow = trackITS->GetNumClust();
+ if(trackITS->GetNumClust()) Chi2Now /= (Double_t )trackITS->GetNumClust();
+ Chi2Ref = reference->GetChi2();
+ Float_t NumClustRef = reference->GetNumClust();
+ if(reference->GetNumClust()) Chi2Ref /= (Double_t )reference->GetNumClust();
+ //cout<<" Chi2Now and Chi2Ref = "<<Chi2Now<<" "<<Chi2Ref<<"\n";
+ if( NumClustNow > NumClustRef ) {*reference = *trackITS;}
+ if((NumClustNow == NumClustRef )&& (Chi2Now < Chi2Ref)) {*reference = *trackITS;}
+ continue;
+ }
+
+ Float_t NumClustNow = trackITS->GetNumClust();
+ if(NumClustNow) {
+ Chi2Now = trackITS->GetChi2();
+ Chi2Now/=NumClustNow;
+ //cout<<" Chi2Now = "<<Chi2Now<<"\n";
+
+ // if(Ptref > 0.6 && Chi2Now > 20.) continue;
+ if(Ptref > 0.6 && Chi2Now > 30.) continue;
+ if((Ptref <= 0.6 && Ptref>0.2)&& Chi2Now > 15.) continue;
+ // if(Chi2Now>5.) continue;
+ //if(Chi2Now>15.) continue;
+ // if(Ptref <= 0.2 && Chi2Now > 10.) continue;
+ if(Ptref <= 0.2 && Chi2Now > 8.) continue;
+ }
+
+ Int_t layerInit = (*trackITS).GetLayer();
+ Int_t layernew = layerInit - 2; // -1 for new layer, -1 for matrix index
+
+ Int_t NLadder[]= {20, 40, 14, 22, 34, 38};
+ Int_t NDetector[]= {4, 4, 5, 8, 23, 26};
+
+ TList listoftrack;
+ Int_t ladp, ladm, detp,detm,ladinters,detinters;
+ Int_t layerfin=layerInit-1;
+ Double_t Rfin=Rlayer[layerfin-1];
+ // cout<<"Prima di intersection \n";
+ Int_t outinters=NewIntersection(*trackITS, Rfin, layerfin, ladinters, detinters);
+
+ // cout<<" outinters = "<<outinters<<"\n";
+ // cout<<" Layer ladder detector intersection ="<<layerfin<<" "<<ladinters<<" "<<detinters<<"\n";
+ // cout << " phiinters zinters = "<<(*trackITS)(0) << " "<<(*trackITS)(1)<<"\n"; getchar();
+
+ if(outinters==-1) continue;
+
+ Int_t flaghit=0;
+
+ if(outinters==0){
+ TVector Touclad(9), Toucdet(9);
+ Int_t lycur=layerfin;
+ ladp=ladinters+1;
+ ladm=ladinters-1;
+ if(ladm <= 0) ladm=NLadder[layerfin-1];
+ if(ladp > NLadder[layerfin-1]) ladp=1;
+ detp=detinters+1;
+ detm=detinters-1;
+ Int_t idetot=1;
+ Touclad(0)=ladinters; Touclad(1)=ladm; Touclad(2)=ladp;
+ Touclad(3)=ladinters; Touclad(4)=ladm; Touclad(5)=ladp;
+ Touclad(6)=ladinters; Touclad(7)=ladm; Touclad(8)=ladp;
+ Toucdet(0)=detinters; Toucdet(1)=detinters; Toucdet(2)=detinters;
+ if(detm > 0 && detp <= NDetector[layerfin-1]) {
+ idetot=9;
+ Toucdet(3)=detm; Toucdet(4)=detm; Toucdet(5)=detm;
+ Toucdet(6)=detp; Toucdet(7)=detp; Toucdet(8)=detp;
+ }
+
+ if(detm > 0 && detp > NDetector[layerfin-1]) {
+ idetot=6;
+ Toucdet(3)=detm; Toucdet(4)=detm; Toucdet(5)=detm;
+ }
+
+ if(detm <= 0 && detp <= NDetector[layerfin-1]) {
+ idetot=6;
+ Toucdet(3)=detp; Toucdet(4)=detp; Toucdet(5)=detp;
+ }
+
+ for (Int_t iriv=0; iriv<idetot; iriv++) { //for on detectors
+ AliITSgeom *g1 = aliITS->GetITSgeom();
+ TVector CTF(9);
+ g1->GetCenterThetaPhi(layerInit-1,(Int_t)Touclad(iriv),(Int_t)Toucdet(iriv),CTF);
+
+ // cout<<" layer, ladder, det, xo, yo, zo = "<<layerInit-1<<" "<<(Int_t)Touclad(iriv)<<
+ // " "<<(Int_t)Toucdet(iriv)<<" "<<CTF(0)<<" "<<CTF(1)<<" "<<CTF(2)<< " "<<CTF(6)<<"\n"; getchar();
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////
+
+ /*** Rec points sorted by module *****/
+ /**************************************/
+
+ Int_t index;
+ AliITSRecPoint *recp;
+ AliITSgeom *geom = aliITS->GetITSgeom();
+ index = geom->GetModuleIndex(lycur,Touclad(iriv),Toucdet(iriv));
+ Int_t lay,lad,det;
+ geom->GetModuleId(index,lay,lad,det);
+ aliITS->ResetRecPoints();
+ gAlice->TreeR()->GetEvent(index+1); //first entry in TreeR is empty
+
+ Int_t npoints=rpoints->GetEntries();
+ Int_t *indlist=new Int_t[npoints+1];
+ Int_t counter=0;
+ Int_t ind;
+ for (ind=0; ind<=npoints; ind++) {
+ indlist[ind]=-1;
+ if (*(vettid[index]+ind)==0) {
+ indlist[counter]=ind;
+ counter++;
+ }
+ }
+
+ ind=-1;
+ for(;;) {
+ ind++;
+ if(indlist[ind] < 0) recp=0;
+ else recp = (AliITSRecPoint*)rpoints->UncheckedAt(indlist[ind]);
+
+ if((!recp) ) break;
+ TVector cluster(3),vecclust(9);
+ vecclust(6)=vecclust(7)=vecclust(8)=-1.;
+ Double_t sigma[2];
+ // set veclust in global
+ Float_t global[3], local[3];
+ local[0]=recp->GetX();
+ local[1]=0.;
+ local[2]= recp->GetZ();
+ AliITSgeom *g1 = aliITS->GetITSgeom();
+ Int_t play = lycur;
+ Int_t plad = TMath::Nint(Touclad(iriv));
+ Int_t pdet = TMath::Nint(Toucdet(iriv));
+ g1->LtoG(play,plad,pdet,local,global);
+
+ vecclust(0)=global[0];
+ vecclust(1)=global[1];
+ vecclust(2)=global[2];
+ vecclust(3) = (float)recp->fTracks[0];
+ vecclust(4) = (float)indlist[ind];
+ vecclust(5) = (float)index;
+ vecclust(6) = (float)recp->fTracks[0];
+ vecclust(7) = (float)recp->fTracks[1];
+ vecclust(8) = (float)recp->fTracks[2];
+
+ sigma[0] = (Double_t) recp->GetSigmaX2();
+ sigma[1] = (Double_t) recp->GetSigmaZ2();
+
+ //now we are in r,phi,z in global
+ cluster(0) = TMath::Sqrt(vecclust(0)*vecclust(0)+vecclust(1)*vecclust(1));//r hit
+ cluster(1) = PhiDef(vecclust(0),vecclust(1)); // phi hit
+ cluster(2) = vecclust(2); // z hit
+ // cout<<" cluster(1) prima = "<<cluster(1)<<"\n";
+ //cluster(1)= cluster(1)-trackITS->Getalphaprov(); //provvisorio;
+ //if(cluster(1)<0.) cluster(1)+=2.*TMath::Pi(); //provvisorio
+ //cout<<" cluster(1) dopo = "<<cluster(1)<< " alphaprov = "<<trackITS->Getalphaprov()<<"\n";
+ Float_t sigmatotphi, sigmatotz;
+
+ //Float_t epsphi=3.2, epsz=3.;
+ Float_t epsphi=3.2, epsz=3.0;
+ //if(Ptref<0.2) {epsphi=3.; epsz=3.;}
+
+ Double_t Rtrack=(*trackITS).Getrtrack();
+ Double_t sigmaphi=sigma[0]/(Rtrack*Rtrack);
+ sigmatotphi=epsphi*TMath::Sqrt(sigmaphi + (*trackITS).GetSigmaphi());
+ sigmatotz=epsz*TMath::Sqrt(sigma[1] + (*trackITS).GetSigmaZ());
+ //cout<<"cluster e sigmatotphi e track = "<<cluster(0)<<" "<<cluster(1)<<" "<<sigmatotphi<<" "<<vecclust(3)<<"\n";
+ //if(vecclust(3)==481) getchar();
+ if(cluster(1)<6. && (*trackITS).Getphi()>6.) cluster(1)=cluster(1)+(2.*TMath::Pi());
+ if(cluster(1)>6. && (*trackITS).Getphi()<6.) cluster(1)=cluster(1)-(2.*TMath::Pi());
+ if(TMath::Abs(cluster(1)-(*trackITS).Getphi()) > sigmatotphi) continue;
+ // cout<<" supero sigmaphi \n";
+ AliITStrack *newTrack = new AliITStrack((*trackITS));
+ (*newTrack).SetLayer((*trackITS).GetLayer()-1);
+
+ if (TMath::Abs(Rtrack-cluster(0))/Rtrack>1e-6)
+ (*newTrack).Correct(Double_t(cluster(0)));
+ //cout<<" cluster(2) e (*newTrack).GetZ() = "<<cluster(2)<<" "<< (*newTrack).GetZ()<<"\n";
+ if(TMath::Abs(cluster(2)-(*newTrack).GetZ()) > sigmatotz){
+ delete newTrack;
+ continue;}
+
+ if(iriv == 0) flaghit=1;
+
+ (*newTrack).AddMS(); // add the multiple scattering matrix to the covariance matrix
+ (*newTrack).AddEL(1.,0);
+ Double_t sigmanew[2];
+ sigmanew[0]= sigmaphi;
+ sigmanew[1]=sigma[1];
+ //cout<<" Chiamo Kalman \n"; getchar();
+ if(flagvert)
+ KalmanFilterVert(newTrack,cluster,sigmanew);
+ else
+ KalmanFilter(newTrack,cluster,sigmanew);
+
+
+ (*newTrack).PutCluster(layernew, vecclust);
+ newTrack->AddClustInTrack();
+
+ listoftrack.AddLast(newTrack);
+
+ } // end of for(;;) on rec points
+
+ delete [] indlist;
+
+ } // end of for on detectors
+ }//end if(outinters==0)
+
+ if(flaghit==0 || outinters==-2) {
+ AliITStrack *newTrack = new AliITStrack(*trackITS);
+ (*newTrack).SetLayer((*trackITS).GetLayer()-1);
+ (*newTrack).AddMS(); // add the multiple scattering matrix to the covariance matrix
+ (*newTrack).AddEL(1.,0);
+
+ listoftrack.AddLast(newTrack);
+ }
+
+
+ //gObjectTable->Print(); // stampa memoria
+
+ AliITStracking(&listoftrack, reference, aliITS, rpoints,Ptref,vettid,flagvert);
+ listoftrack.Delete();
+ } // end of for on tracks
+
+ //gObjectTable->Print(); // stampa memoria
+
+}
+
+
+Int_t AliITStracking::NewIntersection(AliITStrack &track, Double_t rk,Int_t layer, Int_t &ladder, Int_t &detector) {
+//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
+// Found the intersection and the detector
+
+ if(track.DoNotCross(rk)){ /*cout<< " Do not cross \n";*/ return -1;}
+ track.Propagation(rk);
+ Double_t zinters=track.GetZ();
+ Double_t phinters=track.Getphi();
+ //phinters = phinters+track.Getalphaprov(); //provvisorio
+ //if(phinters>2.*3.14) phinters=phinters-2.*3.14; //provvisorio
+ //cout<<"zinters = "<<zinters<<" phinters = "<<phinters<<"\n";
+
+ ////////////////////////////////// limits for Geometry 5 /////////////////////////////
+
+ Int_t NLadder[]= {20, 40, 14, 22, 34, 38};
+ Int_t NDetector[]= {4, 4, 5, 8, 23, 26};
+
+ Float_t Detx[]= {0.64, 0.64, 3.509, 3.509, 3.65, 3.65 };
+ Float_t Detz[]= {4.19, 4.19, 3.75 , 3.75 , 2 , 2 };
+ ////////////////////////////////////////////////////////////////////////////////////////////////
+
+ TVector det(9);
+ TVector ListDet(2);
+ TVector DistZCenter(2);
+ AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
+ Int_t iz=0;
+ Double_t epsz=1.2;
+ Double_t epszpixel=0.05;
+
+ for(Int_t iD = 1; iD<= NDetector[layer-1]; iD++) {
+ g1->GetCenterThetaPhi(layer,1,iD,det);
+ Double_t zmin=det(2)-Detz[layer-1];
+ if(iD==1) zmin=det(2)-(Detz[layer-1])*epsz;
+ Double_t zmax=det(2)+Detz[layer-1];
+ if(iD==NDetector[layer-1]) zmax=det(2)+(Detz[layer-1])*epsz;
+ if(layer == 3 || layer==2) zmin=zmin-epszpixel; zmax=zmax+epszpixel;
+ //cout<<"zmin zinters zmax det(2)= "<<zmin<<" "<<zinters<<" "<<zmax<<" "<<det(2)<<"\n";
+ if(zinters > zmin && zinters <= zmax) {
+ if(iz>1) {cout<< " Errore su iz in NewIntersection \n"; getchar();}
+ else {
+ ListDet(iz)= iD; DistZCenter(iz)=TMath::Abs(zinters-det(2)); iz++;
+ }
+ }
+ }
+
+ if(iz==0) {/* cout<< " No detector along Z \n";*/ return -2;}
+ detector=Int_t (ListDet(0));
+ if(iz>1 && (DistZCenter(0)>DistZCenter(1))) detector=Int_t (ListDet(1));
+
+ AliITSgeom *g2 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
+ Float_t global[3];
+ Float_t local[3];
+ TVector ListLad(2);
+ TVector DistphiCenter(2);
+ Int_t ip=0;
+ Double_t pigre=TMath::Pi();
+
+ for(Int_t iLd = 1; iLd<= NLadder[layer-1]; iLd++) {
+ g1->GetCenterThetaPhi(layer,iLd,detector,det);
+ Double_t phidet=PhiDef(Double_t(det(0)),Double_t(det(1)));
+
+
+ // cout<<" layer phidet e det(6) = "<< layer<<" "<<phidet<<" "<<det(6)<<"\n"; getchar();
+ Double_t xmin,ymin,xmax,ymax;
+ Double_t phiconfr;
+ local[1]=local[2]=0.;
+ local[0]= -(Detx[layer-1]);
+ if(layer==1) local[0]= (Detx[layer-1]); //take into account different reference system
+ g2->LtoG(layer,iLd,detector,local,global);
+ xmax=global[0]; ymax=global[1];
+ local[0]= (Detx[layer-1]);
+ if(layer==1) local[0]= -(Detx[layer-1]); //take into account different reference system
+ g2->LtoG(layer,iLd,detector,local,global);
+ xmin=global[0]; ymin=global[1];
+ Double_t phimin=PhiDef(xmin,ymin);
+ Double_t phimax=PhiDef(xmax,ymax);
+ //cout<<" xmin ymin = "<<xmin<<" "<<ymin<<"\n";
+ // cout<<" xmax ymax = "<<xmax<<" "<<ymax<<"\n";
+ // cout<<" iLd phimin phimax ="<<iLd<<" "<<phimin<<" "<<phimax<<"\n"; getchar();
+ if(phimin>phimax ){
+ if(phimin <5.5) {cout<<" Error in NewIntersection for phi \n"; getchar();}
+ phimin=phimin-(2.*pigre);
+ if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre);
+ if(phidet>(1.5*pigre)) phidet=phidet-(2.*pigre);
+ }
+ else phiconfr=phinters;
+ if(phiconfr>phimin && phiconfr<= phimax) {
+ if(ip>1) {
+ cout<< " Errore su ip in NewIntersection \n"; getchar();
+ }
+ else {
+ ListLad(ip)= iLd; DistphiCenter(ip)=TMath::Abs(phiconfr-phidet); ip++;
+ }
+ }
+ }
+ if(ip==0) { cout<< " No detector along phi \n"; getchar();}
+ ladder=Int_t (ListLad(0));
+ if(ip>1 && (DistphiCenter(0)>DistphiCenter(1))) ladder=Int_t (ListLad(1));
+
+ return 0;
+}
+
+
+Double_t AliITStracking::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;
+}
+
+
+void AliITStracking::KalmanFilter(AliITStrack *newTrack,TVector &cluster,Double_t sigma[2]){
+//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
+// Kalman filter without vertex constraint
+
+ TMatrix H(2,5); H.UnitMatrix();
+ TMatrix Ht(TMatrix::kTransposed, H);
+
+
+ ////////////////////////////// Evaluation of the measurement vector /////////////////////////////////////
+
+ TVector m(2);
+ Double_t rk,phik,zk;
+ rk=cluster(0); phik=cluster(1); zk=cluster(2);
+ m(0)=phik; m(1)=zk;
+ // cout<<" r and m = "<<rk<<" "<<m(0)<<" "<<m(1)<<"\n";
+ ///////////////////////////////////// Evaluation of the error matrix V ///////////////////////////////
+
+ TMatrix V(2,2);
+ V(0,1)=0.; V(1,0)=0.;
+ V(0,0)=sigma[0];
+ V(1,1)=sigma[1];
+ ///////////////////////////////////////////////////////////////////////////////////////////
+
+ TMatrix C=newTrack->GetCMatrix();
+ TMatrix tmp(H,TMatrix::kMult,C);
+ TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
+
+ R.Invert();
+ // cout<<" R prima = \n";
+ // R.Print(); getchar();
+
+ TMatrix K(C,TMatrix::kMult,Ht); K*=R;
+
+ TVector x=newTrack->GetVector();
+
+ TVector savex=x;
+ x*=H; x-=m;
+ x*=-1.; x*=K; x+=savex;
+ TMatrix saveC=C;
+ C.Mult(K,tmp); C-=saveC; C*=-1;
+ newTrack->GetVector()=x;
+ newTrack->GetCMatrix()=C;
+
+ TVector res= newTrack->GetVector();
+ //cout<<" res = "<<res(0)<<" "<<res(1)<<" "<<res(2)<<" "<<res(3)<<" "<<res(4)<<"\n";
+ res*=H; res-=m; res*=-1.;
+ TMatrix Cn=newTrack->GetCMatrix();
+ TMatrix tmpn(H,TMatrix::kMult,Cn);
+ TMatrix Rn(tmpn,TMatrix::kMult,Ht); Rn-=V; Rn*=-1.;
+
+ Rn.Invert();
+ TVector r=res; res*=Rn;
+ //cout<<" R dopo = \n";
+ //Rn.Print(); getchar();
+ Double_t chi2= r*res;
+ //cout<<"chi2 ="<<chi2<<"\n"; getchar();
+
+ newTrack->SetChi2(newTrack->GetChi2()+chi2);
+
+}
+
+
+void AliITStracking::KalmanFilterVert(AliITStrack *newTrack,TVector &cluster,Double_t sigma[2]){
+//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
+// Kalman filter with vertex constraint
+ TMatrix H(4,5); H.UnitMatrix();
+ TMatrix Ht(TMatrix::kTransposed, H);
+
+ ////////////////////////////// Evaluation of the measurement vector /////////////////////////////////////
+
+ TVector m(4);
+ Double_t rk,phik,zk;
+ rk=cluster(0); phik=cluster(1); zk=cluster(2);
+ m(0)=phik; m(1)=zk;
+
+ Double_t CC=(*newTrack).GetC();
+ Double_t Zv=(*newTrack).GetZv();
+ Double_t Dv=(*newTrack).GetDv();
+ // cout<<" Dv e Zv = "<<Dv<<" "<<Zv<<"\n";
+ Double_t Cy=CC/2.;
+ Double_t tgl= (zk-Zv)*Cy/TMath::ASin(Cy*rk);
+ m(2)=Dv; m(3)=tgl;
+ // cout<<" m = \n";
+ // cout<<m(0)<<" "<<m(1)<<" "<<m(2)<<" "<<m(3)<<"\n";
+ ///////////////////////////////////// Evaluation of the error matrix V ///////////////////////////////
+
+ TMatrix V(4,4);
+ V(0,1)=0.; V(1,0)=0.;
+ Int_t Layer=newTrack->GetLayer();
+ V(0,0)=sigma[0];
+ V(1,1)=sigma[1];
+ V(1,3)=sigma[1]/rk; V(3,1)=V(1,3);
+ Double_t sigmaDv=newTrack->GetsigmaDv();
+ V(2,2)=sigmaDv*sigmaDv + newTrack->Getd2(Layer-1);
+ V(2,3)=newTrack->Getdtgl(Layer-1); V(3,2)=V(2,3);
+ Double_t sigmaZv=newTrack->GetsigmaZv();
+ V(3,3)=(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(Layer-1);
+ ///////////////////////////////////////////////////////////////////////////////////////////
+
+ //cout<<" d2 tgl2 dtgl = "<<newTrack->Getd2(Layer-1)<<" "<<newTrack->Gettgl2(Layer-1)<<" "<<newTrack->Getdtgl(Layer-1)<<"\n";
+ // cout<<" V = \n";
+ // V.Print(); getchar();
+
+ TMatrix C=newTrack->GetCMatrix();
+ TMatrix tmp(H,TMatrix::kMult,C);
+ TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
+
+ R.Invert();
+ TMatrix K(C,TMatrix::kMult,Ht); K*=R;
+ TVector x=newTrack->GetVector();
+ TVector savex=x;
+ x*=H; x-=m;
+ x*=-1; x*=K; x+=savex;
+ TMatrix saveC=C;
+ C.Mult(K,tmp); C-=saveC; C*=-1;
+ newTrack->GetVector()=x;
+ newTrack->GetCMatrix()=C;
+ TVector res= newTrack->GetVector();
+ //cout<<" res = "<<res(0)<<" "<<res(1)<<" "<<res(2)<<" "<<res(3)<<" "<<res(4)<<"\n";
+ res*=H; res-=m; res*=-1.;
+ TMatrix Cn=newTrack->GetCMatrix();
+ TMatrix tmpn(H,TMatrix::kMult,Cn);
+ TMatrix Rn(tmpn,TMatrix::kMult,Ht); Rn-=V; Rn*=-1.;
+
+ Rn.Invert();
+ TVector r=res; res*=Rn;
+ Double_t chi2= r*res;
+
+ newTrack->SetChi2(newTrack->GetChi2()+chi2);
+
+}
+
--- /dev/null
+#ifndef ALIITSTRACKING_H
+#define ALIITSTRACKING_H
+
+#include <TObject.h>
+#include <TList.h>
+
+class TObjArray;
+class TVector;
+class TMatrix;
+class AliITStrack;
+class AliITS;
+
+
+class AliITStracking : public TObject {
+
+//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
+
+ Double_t Rlayer[6];
+
+public:
+
+ AliITStracking() {;}
+
+ AliITStracking(TList *trackITSlist,AliITStrack *reference,AliITS *obj,TObjArray *fpoints,
+ Double_t Ptref, Int_t **vettid, Bool_t flagvert );
+
+ Int_t NewIntersection(AliITStrack &track, Double_t rk,Int_t layer, Int_t &ladder, Int_t &detector );
+ Double_t PhiDef(Double_t x, Double_t y);
+
+ void KalmanFilter(AliITStrack *newtrack, TVector &cluster, Double_t sigma[2]);
+ void KalmanFilterVert(AliITStrack *newtrack, TVector &cluster, Double_t sigma[2]);
+
+ ClassDef(AliITStracking,1)
+};
+
+#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff