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