ITS/AliITStrack.cxx

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