ITS/AliITSTrackV1.cxx

   1 //   ITS Track Class
   2 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
   3 // It contain all the usefull information for the track and the method to calculate, modify or extract them
   4 // The track is mainly caracterized by the state vector of elements (fX0, fX1, fX2, fX3, fX4) and the
   5 // corresponding covariance matrix of elements (C00, C10, ..... C44) that is triangular
   6 //
   7 #include <iostream.h>
   8 #include <TMath.h>
   9 #include <TVector.h>
  10 #include <TObjArray.h>
  11 #include "AliRun.h"
  12 #include "AliITSRad.h"
  13 #include "AliITSTrackV1.h"
  14 #include "AliGenerator.h"
  15 #include "AliMagF.h"
  16
  17
  18 ClassImp(AliITSTrackV1)
  19
  20 AliITSTrackV1::AliITSTrackV1() {
  21 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
  22 // default constructor
  23
  24   fTPCtrack = 0;
  25   fC00=fC10=fC11=fC20=fC21=fC22=fC30=fC31=fC32=fC33=fC40=fC41=fC42=fC43=fC44=0.;
  26   flistCluster = new TObjArray;
  27   fNumClustInTrack =0;
  28   fChi2=-1;
  29   flabel =0;
  30   fVertex.ResizeTo(3);
  31   fErrorVertex.ResizeTo(3);
  32   fLayer = -1;
  33   fClusterInTrack = new TMatrix(6,9);
  34   Int_t i,j;
  35   //for(i=0; i<6; i++) (*fClusterInTrack)(i,6)=(*fClusterInTrack)(i,7)=
  36   //                         (*fClusterInTrack)(i,8)=-1.;
  37   for(i=0; i<6; i++){
  38   for(j=4; j<9; j++) (*fClusterInTrack)(i,j)=-1.;   //modificata angela
  39   }
  40   frtrack=0.;
  41   fnoclust=0;
  42   fd2.ResizeTo(6);
  43   ftgl2.ResizeTo(6);
  44   fdtgl.ResizeTo(6);
  45
  46 //////////////////////////////////////// gets magnetic field factor ////////////////////////////////
  47
  48   AliMagF * fieldPointer = gAlice->Field();
  49   fFieldFactor = (Double_t)fieldPointer->Factor();
  50   //cout<< " field factor = "<<fFieldFactor<<"\n"; getchar();
  51
  52 /////////////////////////////////////////////////////////////////////////////////////////////////////////
  53
  54 }
  55
  56
  57
  58 AliITSTrackV1::AliITSTrackV1(const AliITSTrackV1 &cobj) {
  59 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
  60 // copy constructor
  61
  62   fClusterInTrack = new TMatrix(6,9);
  63   Int_t i,j;
  64   //for(i=0; i<6; i++) (*fClusterInTrack)(i,6)=(*fClusterInTrack)(i,7)=
  65   //                         (*fClusterInTrack)(i,8)=-1.;
  66   for(i=0; i<6; i++){
  67   for(j=4; j<9; j++) (*fClusterInTrack)(i,j)=-1.;   //modificata angela
  68   }
  69   flistCluster = new TObjArray;
  70   fVertex.ResizeTo(3);
  71   fErrorVertex.ResizeTo(3);
  72   fVertex = cobj.fVertex;
  73   fErrorVertex = cobj.fErrorVertex;
  74   flabel = cobj.flabel;
  75   fLayer=cobj.fLayer;
  76   fTPCtrack = cobj.fTPCtrack;
  77   fNumClustInTrack = cobj.fNumClustInTrack;
  78   fChi2= cobj.fChi2;
  79   fX0=cobj.fX0; fX1=cobj.fX1; fX2=cobj.fX2; fX3=cobj.fX3; fX4=cobj.fX4;
  80   frtrack=cobj.frtrack;
  81   fDv=cobj.fDv;
  82   fZv=cobj.fZv;
  83   fsigmaDv=cobj.fsigmaDv;
  84   fsigmaZv=cobj.fsigmaZv;
  85   fd2.ResizeTo(6);
  86   ftgl2.ResizeTo(6);
  87   fdtgl.ResizeTo(6);
  88   fd2=cobj.fd2;
  89   ftgl2=cobj.ftgl2;
  90   fdtgl=cobj.fdtgl;
  91   fnoclust=cobj.fnoclust;
  92
  93   fC00=cobj.fC00; fC10=cobj.fC10; fC11=cobj.fC11; fC20=cobj.fC20; fC21=cobj.fC21;
  94   fC22=cobj.fC22; fC30=cobj.fC30; fC31=cobj.fC31; fC32=cobj.fC32; fC33=cobj.fC33;
  95   fC40=cobj.fC40; fC41=cobj.fC41; fC42=cobj.fC42; fC43=cobj.fC43; fC44=cobj.fC44;
  96
  97   *fClusterInTrack = *cobj.fClusterInTrack;
  98
  99   fFieldFactor=cobj.fFieldFactor;
 100
 101   for(i=0; i<cobj.flistCluster->GetSize(); i++)
 102     flistCluster->AddLast(cobj.flistCluster->At(i));
 103
 104 }
 105
 106 AliITSTrackV1::AliITSTrackV1(AliTPCtrack &obj)
 107 {
 108 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 109 // special constructor to convert a TPC track into an ITS track
 110
 111 //////////////////////////////////////// gets magnetic field factor ////////////////////////////////
 112
 113   AliMagF * fieldPointer = gAlice->Field();
 114   fFieldFactor = (Double_t)fieldPointer->Factor();
 115   //cout<< " field factor = "<<fFieldFactor<<"\n"; getchar();
 116
 117 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 118
 119
 120   fTPCtrack = &obj;
 121   fVertex.ResizeTo(3);
 122   fErrorVertex.ResizeTo(3);
 123
 124   fd2.ResizeTo(6);
 125   ftgl2.ResizeTo(6);
 126   fdtgl.ResizeTo(6);
 127   AliGenerator *gener = gAlice->Generator();
 128   Float_t vxg,vyg,vzg;
 129   gener->GetOrigin(vxg,vyg,vzg);
 130
 131
 132   fVertex(0)=(Double_t)vxg;
 133   fVertex(1)=(Double_t)vyg;
 134   fVertex(2)=(Double_t)vzg;
 135
 136   fLayer = 7;
 137   //fmCovariance = new TMatrix(5,5);
 138   fClusterInTrack = new TMatrix(6,9);
 139
 140   Int_t i,j;
 141   //for(i=0; i<6; i++) (*fClusterInTrack)(i,6)=(*fClusterInTrack)(i,7)=
 142   //                         (*fClusterInTrack)(i,8)=-1.;
 143   for(i=0; i<6; i++){
 144   for(j=4; j<9; j++) (*fClusterInTrack)(i,j)=-1.;   //modificata angela
 145   }
 146   flistCluster = new TObjArray;
 147   fNumClustInTrack = 0;
 148   fnoclust=0;
 149   LmTPC();
 150
 151 }
 152
 153
 154 AliITSTrackV1::~AliITSTrackV1() {
 155 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 156 //destructor
 157
 158   if(flistCluster) {
 159     flistCluster->Delete();
 160     delete flistCluster;
 161   }
 162   if(fClusterInTrack) delete fClusterInTrack;
 163
 164 }
 165
 166 void AliITSTrackV1::PutCElements(Double_t C00, Double_t C10, Double_t C11, Double_t C20, Double_t C21,
 167 Double_t C22, Double_t C30, Double_t C31, Double_t C32, Double_t C33, Double_t C40,
 168 Double_t C41, Double_t C42, Double_t C43, Double_t C44){
 169 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 170 // puts elements of covariance matrix
 171
 172   fC00=C00; fC10=C10; fC11=C11; fC20=C20; fC21=C21; fC22=C22; fC30=C30; fC31=C31; fC32=C32; fC33=C33;
 173   fC40=C40; fC41=C41; fC42=C42; fC43=C43; fC44=C44;
 174 }
 175
 176 void AliITSTrackV1::GetCElements(Double_t &C00, Double_t &C10, Double_t &C11, Double_t &C20, Double_t &C21,
 177 Double_t &C22, Double_t &C30, Double_t &C31, Double_t &C32, Double_t &C33, Double_t &C40,
 178 Double_t &C41, Double_t &C42, Double_t &C43, Double_t &C44) const {
 179 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 180 // gets elements of covariance matrix
 181
 182
 183   C00=fC00; C10=fC10; C11=fC11; C20=fC20; C21=fC21; C22=fC22; C30=fC30; C31=fC31; C32=fC32; C33=fC33;
 184   C40=fC40; C41=fC41; C42=fC42; C43=fC43; C44=fC44;
 185
 186 }
 187
 188 void AliITSTrackV1::GetXElements(Double_t &X0, Double_t &X1, Double_t &X2, Double_t &X3, Double_t &X4) const {
 189 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 190 // gets the elements of the state vector
 191   X0=fX0; X1=fX1; X2=fX2; X3=fX3; X4=fX4;
 192 }
 193
 194 void AliITSTrackV1::PutXElements(Double_t X0, Double_t X1, Double_t X2, Double_t X3, Double_t X4){
 195 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 196 // puts the elements of the state vector
 197   fX0=X0; fX1=X1; fX2=X2; fX3=X3; fX4=X4;
 198 }
 199
 200 void AliITSTrackV1::LmTPC() {
 201 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 202 // Transform the TPC state vector from TPC-local to master and build a new state vector ITS-type
 203 // The covariance matrix is also modified accordingly
 204
 205
 206   Double_t alpha = fTPCtrack->GetAlpha();
 207
 208   //printf("LmTPC: alpha %f\n",alpha);
 209
 210   Double_t yTPC   = fTPCtrack->GetY();
 211   Double_t zTPC   = fTPCtrack->GetZ();
 212   Double_t cTPC   = fTPCtrack->GetC();
 213   Double_t etaTPC = fTPCtrack->GetEta();
 214   Double_t tglTPC = fTPCtrack->GetTgl();
 215
 216
 217   Double_t xm, ym, zm;
 218   Double_t sina = TMath::Sin(alpha);
 219   Double_t cosa = TMath::Cos(alpha);
 220   Double_t xl= fTPCtrack->GetX();
 221   xm = xl * cosa - yTPC*sina;
 222   ym = xl * sina + yTPC*cosa;
 223   zm = zTPC;
 224   //cout<<" xl e alpha = "<<xl<<" "<<alpha<<"\n"; getchar();
 225
 226   Double_t x0m,y0m;
 227
 228   ///////////////////////////////////// determine yo //////////////////////////////////////////////////
 229
 230   Double_t vxl=fVertex(0)*cosa+fVertex(1)*sina;
 231   Double_t vyl= -fVertex(0)*sina+fVertex(1)*cosa;
 232   Double_t xo,yo, signy;
 233   Double_t r = 1./cTPC;
 234   xo =  etaTPC / cTPC;
 235   //fxoTPC=xo;
 236   Double_t yo1, yo2, diffsq1, diffsq2;
 237   yo1 = yTPC +  TMath::Sqrt(r*r - (xl-xo)*(xl-xo));
 238   yo2 = yTPC -  TMath::Sqrt(r*r - (xl-xo)*(xl-xo));
 239   diffsq1=TMath::Abs((yo1-vyl)*(yo1-vyl)+(xo-vxl)*(xo-vxl)-r*r);
 240   diffsq2=TMath::Abs((yo2-vyl)*(yo2- vyl)+(xo-vxl)*(xo-vxl)-r*r);
 241   if(diffsq1<diffsq2) {yo=yo1; signy=1.;} else {yo=yo2; signy=-1.;};
 242
 243   ////////////////////////////////////////////////////////////////////////////////////////////////////
 244
 245   x0m = xo * cosa - yo * sina;
 246   y0m = xo * sina + yo * cosa;
 247
 248   frtrack=TMath::Sqrt(xm*xm+ym*ym);
 249   /*
 250   Double_t pigre=TMath::Pi();
 251   Double_t phi=0.0;
 252   if(ym == 0. || xm == 0.) {
 253     if(ym == 0. && xm == 0.) {cout << "  Error in AliITSTrackV1::LmTPC x=0 and y=0 \n"; getchar();}
 254     if(ym ==0. && xm>0.) phi=0.;
 255     if(ym==0. && xm<0.) phi=pigre;
 256     if(xm==0 && ym>0.) phi=pigre/2.;
 257     if(xm==0 && ym<0.) phi=1.5*pigre;
 258   }
 259   else {
 260     if (xm>0. && ym>0.) phi=TMath::ATan(ym/xm);
 261     if (xm<0. && ym>0.) phi=pigre+TMath::ATan(ym/xm);
 262     if (xm<0. && ym<0.) phi=pigre+TMath::ATan(ym/xm);
 263     if (xm>0. && ym<0.) phi=(2*pigre)+TMath::ATan(ym/xm);
 264   };
 265   if(phi<0. || phi>(2*pigre)) {cout<<"attention error on phi in  AliITSTrackV1:LmTPC \n"; getchar();}
 266    */
 267
 268     Double_t phi=TMath::ATan2(ym,xm);  if(phi<0) phi=2.*TMath::Pi()+phi;   //nuova def phi
 269
 270   fX0=phi;
 271   fX1=zm;
 272   fX3=tglTPC;
 273   fX4=cTPC;
 274
 275
 276   Double_t dd=TMath::Sqrt((x0m-fVertex(0))*(x0m-fVertex(0))+(y0m-fVertex(1))*(y0m-fVertex(1)));
 277   Double_t signdd;
 278   if (r>0) signdd=1.; else signdd=-1.;
 279   fX2=signdd*dd-r;
 280   //cout<<" fvertex = "<<fVertex(0)<<" "<<fVertex(1)<<" "<<fVertex(2)<<"\n";
 281
 282   Double_t cov[15];
 283   fTPCtrack->GetCovariance(cov);
 284
 285   Double_t  dfidy, dDdy, dDdC, dDdeta;
 286
 287   dfidy=(xm*cosa+ym*sina)/(frtrack*frtrack);
 288   dDdy=signdd*((y0m-fVertex(1))*cosa-(x0m-fVertex(0))*sina)/dd;
 289   Double_t dyodr=signy*(r+(xl-xo)*etaTPC)/TMath::Sqrt(r*r-(xl-xo)*(xl-xo));
 290   Double_t dyomdr=sina*etaTPC+cosa*dyodr;
 291   Double_t dxomdr=cosa*etaTPC-sina*dyodr;
 292   Double_t ddddR=((x0m-fVertex(0))*dxomdr+(y0m-fVertex(1))*dyomdr)/dd;
 293   dDdC=-r*r*(signdd*ddddR-1.);
 294   Double_t dyoldxol=signy*(xl-xo)/TMath::Sqrt(r*r-(xl-xo)*(xl-xo));
 295   Double_t dxomdeta=r*(cosa-sina*dyoldxol);
 296   Double_t dyomdeta=r*(sina+cosa*dyoldxol);
 297   dDdeta=signdd*((x0m-fVertex(0))*dxomdeta+(y0m-fVertex(1))*dyomdeta)/dd;
 298
 299   Double_t f00=dfidy;
 300   Double_t f20=dDdy;
 301   Double_t f22=dDdC;
 302   Double_t f23=dDdeta;
 303
 304   Double_t t00=cov[0]*f00;
 305   Double_t t02=cov[0]*f20+cov[6]*f22+cov[3]*f23;
 306   Double_t t20=cov[6]*f00;
 307   Double_t t22=cov[6]*f20+cov[9]*f22+cov[8]*f23;
 308
 309   fC00=f00*t00;
 310   fC10=cov[1]*f00;
 311   fC11=cov[2];
 312   fC20=f20*t00+f22*t20+f23*cov[3]*f00;
 313   fC21=f20*cov[1]+f22*cov[7]+f23*cov[4];
 314   fC22=f20*t02+f22*t22+f23*(cov[3]*f20+cov[8]*f22+cov[5]*f23);
 315   fC30=cov[10]*f00;
 316   fC31=cov[11];
 317   fC32=cov[10]*f20+cov[13]*f22+cov[12]*f23;
 318   fC33=cov[14];
 319   fC40=t20;
 320   fC41=cov[7];
 321   fC42=t22;
 322   fC43=cov[13];
 323   fC44=cov[9];
 324
 325   //cout<<" C32 e C44 = "<<fC32<<" "<<fC44<<"\n"; getchar();
 326
 327 }
 328
 329
 330 AliITSTrackV1 &AliITSTrackV1::operator=(AliITSTrackV1 obj) {
 331 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 332 // assignement operator
 333
 334   if (flistCluster) {
 335     flistCluster->Delete();
 336     delete flistCluster;
 337   }
 338   delete fClusterInTrack;
 339   fClusterInTrack = new TMatrix(6,9);
 340   flistCluster = new TObjArray;
 341   flabel = obj.flabel;
 342   fTPCtrack = obj.fTPCtrack;
 343   fNumClustInTrack = obj.fNumClustInTrack;
 344   fChi2= obj.fChi2;
 345   fVertex=obj.fVertex;
 346   fErrorVertex=obj.fErrorVertex;
 347   fX0=obj.fX0; fX1=obj.fX1; fX2=obj.fX2; fX3=obj.fX3; fX4=obj.fX4;
 348   fLayer=obj.fLayer;
 349   frtrack=obj.frtrack;
 350   fDv=obj.fDv;
 351   fZv=obj.fZv;
 352   fsigmaDv=obj.fsigmaDv;
 353   fsigmaZv=obj.fsigmaZv;
 354   fd2=obj.fd2;
 355   ftgl2=obj.ftgl2;
 356   fdtgl=obj.fdtgl;
 357   fnoclust=obj.fnoclust;
 358
 359   fC00=obj.fC00; fC10=obj.fC10; fC11=obj.fC11; fC20=obj.fC20; fC21=obj.fC21;
 360   fC22=obj.fC22; fC30=obj.fC30; fC31=obj.fC31; fC32=obj.fC32; fC33=obj.fC33;
 361   fC40=obj.fC40; fC41=obj.fC41; fC42=obj.fC42; fC43=obj.fC43; fC44=obj.fC44;
 362
 363
 364   *fClusterInTrack = *obj.fClusterInTrack;
 365   Int_t i;
 366   for(i=0; i<obj.flistCluster->GetSize(); i++) flistCluster->AddLast(obj.flistCluster->At(i));
 367
 368   return *this;
 369
 370 }
 371
 372 void AliITSTrackV1::PutCluster(Int_t layerc, TVector vecclust) {
 373 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 374 // put information for clusters
 375
 376   (*fClusterInTrack)(layerc,0) = vecclust(0);
 377   (*fClusterInTrack)(layerc,1) = vecclust(1);
 378   (*fClusterInTrack)(layerc,2) = vecclust(2);
 379   (*fClusterInTrack)(layerc,3) = vecclust(3);
 380   (*fClusterInTrack)(layerc,4) = vecclust(4);
 381   (*fClusterInTrack)(layerc,5) = vecclust(5);
 382   (*fClusterInTrack)(layerc,6) = vecclust(6);
 383   (*fClusterInTrack)(layerc,7) = vecclust(7);
 384   (*fClusterInTrack)(layerc,8) = vecclust(8);
 385
 386 }
 387
 388
 389 void AliITSTrackV1::GetClusters() {
 390 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 391 // prints the clusters belonging to the current track
 392
 393   TMatrix a(*fClusterInTrack);
 394   TMatrix b(6,3);
 395   Int_t i;
 396   for(i=0;i<6; i++){
 397    b(i,0)=a(i,6); b(i,1)=a(i,7); b(i,2)=a(i,8);
 398   }
 399   a.Print();
 400  // b.Print();
 401
 402 }
 403
 404
 405 TVector AliITSTrackV1::GetLabTrack(Int_t lay) {
 406 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 407 // gets the label of the track
 408
 409   TVector vecLabel(3);
 410   vecLabel(0)=( (Float_t) (*fClusterInTrack)(lay,6) );
 411   vecLabel(1)=( (Float_t) (*fClusterInTrack)(lay,7) );
 412   vecLabel(2)=( (Float_t) (*fClusterInTrack)(lay,8) );
 413   return vecLabel;
 414 }
 415
 416 void AliITSTrackV1::Search(TVector VecTotLabref, Long_t &labref, Int_t &freq){
 417 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 418 // define label
 419
 420   Int_t vecfreq[18];
 421
 422   Int_t i,j;
 423   for(i=0; i<18; i++) vecfreq[i]=0;
 424
 425   for(i=0; i<18; i++) {
 426     for(j=0; j<18; j++) {
 427      // if(VecTotLabref(i) == 0.) VecTotLabref(i)=-3.;  //commentato il 5-3-2001
 428       if( (VecTotLabref(i)>=0.) && (VecTotLabref(i)==VecTotLabref(j)) ) vecfreq[i]++;
 429     }
 430   }
 431   Int_t imax=-1000;
 432   Long_t  labdefault= (Long_t)1000000.;
 433   freq=0;
 434   for(i=0; i<18; i++) {
 435     if(vecfreq[i]>freq) {freq=vecfreq[i]; imax=i;}
 436   }
 437   if(imax<0) labref=labdefault; else labref=(Long_t) VecTotLabref(imax);
 438 }
 439
 440
 441 void AliITSTrackV1::Propagation(Double_t rk) {
 442 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 443 //Propagation of track
 444   Double_t duepi=2.*TMath::Pi();
 445   Double_t rkm1=frtrack;
 446 //cout<<" rk e rkm1 dentro Propagation   "<<rk<<" "<<rkm1<<"\n";
 447
 448   //
 449   Double_t aAk=ArgA(rk), aAkm1=ArgA(rkm1);
 450   Double_t ak=Arga(rk), akm1=Arga(rkm1);
 451   fX0+=TMath::ASin(aAk)-TMath::ASin(aAkm1);
 452
 453   if(fX0>duepi) fX0-=duepi;
 454   if(fX0<0.) fX0+=duepi;
 455
 456   Double_t tgl=fX3;
 457   Double_t c=fX4;
 458   Double_t d=fX2;
 459   Double_t cy=c/2;
 460   fX1+=tgl/cy*(TMath::ASin(ak)-TMath::ASin(akm1));
 461   frtrack=rk;
 462
 463
 464   Double_t bk=ArgB(rk), bkm1=ArgB(rkm1);
 465   Double_t ck=ArgC(rk), ckm1=ArgC(rkm1);
 466
 467   Double_t f02=ck/TMath::Sqrt(1.-aAk*aAk) - ckm1/TMath::Sqrt(1.-aAkm1*aAkm1);
 468   Double_t f04=bk/TMath::Sqrt(1.-aAk*aAk) - bkm1/TMath::Sqrt(1.-aAkm1*aAkm1);
 469   Double_t f12=tgl*d*(1./rk - 1./rkm1);
 470   Double_t f13=rk - rkm1;
 471
 472
 473   Double_t c00=fC00;
 474   Double_t c10=fC10;
 475   Double_t c11=fC11;
 476   Double_t c20=fC20;
 477   Double_t c21=fC21;
 478   Double_t c22=fC22;
 479   Double_t c30=fC30;
 480   Double_t c31=fC31;
 481   Double_t c32=fC32;
 482   Double_t c33=fC33;
 483   Double_t c40=fC40;
 484   Double_t c41=fC41;
 485   Double_t c42=fC42;
 486   Double_t c43=fC43;
 487   Double_t c44=fC44;
 488
 489   Double_t r10=c10+c21*f02+c41*f04;
 490   Double_t r20=c20+c22*f02+c42*f04;
 491   Double_t r30=c30+c32*f02+c43*f04;
 492   Double_t r40=c40+c42*f02+c44*f04;
 493   Double_t r21=c21+c22*f12+c32*f13;
 494   Double_t r31=c31+c32*f12+c33*f13;
 495   Double_t r41=c41+c42*f12+c43*f13;
 496
 497   fC00=c00+c20*f02+c40*f04+f02*r20+f04*r40;
 498   fC10=r10+f12*r20+f13*r30;
 499   fC11=c11+c21*f12+c31*f13+f12*r21+f13*r31;
 500   fC20=r20;
 501   fC21=r21;
 502   fC30=r30;
 503   fC31=r31;
 504   fC40=r40;
 505   fC41=r41;
 506
 507 }
 508
 509 void AliITSTrackV1::AddEL(AliITSRad *rl, Double_t signdE, Bool_t flagtot, Double_t mass) {
 510 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 511 //  add energy loss
 512
 513   TVector s(6);
 514   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;
 515   //s(0)=0.0026+0.00283*2.; s(1)=0.018*2.; s(2)=0.0094; s(3)=0.0095; s(4)=0.0091; s(5)=0.0087;
 516   //0.00277 is added in the first layer to take into account the energy loss in the beam pipe
 517
 518   //for(int k=0; k<6; k++) cout<<s(k)<<" "; cout<<"\n";
 519   Int_t k;
 520   //for(k=0; k<6; k++) s(k)=s(k)*1.6;    //forint
 521   for(k=0; k<6; k++) s(k)=s(k)*1.7;    //forint
 522
 523   Double_t phi=fX0;
 524
 525   if(phi<0.174 ) s(5)=s(5)+0.012;
 526   if(phi>6.1 ) s(5)=s(5)+0.012; // to take into account rail
 527   if(phi>2.96 && phi<3.31 ) s(5)=s(5)+0.012;
 528
 529   /*
 530   Double_t tgl=fX3;
 531   Double_t theta=((TMath::Pi())/2.)-TMath::ATan(tgl);
 532   //phi*=180./TMath::Pi();
 533   //theta*=180./TMath::Pi();
 534   //Double_t rad90=(TMath::Pi())/2.;
 535   Double_t rad40=(TMath::Pi())*40./180.;
 536   Double_t rad100=(TMath::Pi())*100/180;
 537   Double_t rad360=(TMath::Pi())*2.;
 538   Int_t imax=rl->Getimax();
 539   Int_t jmax=rl->Getjmax();
 540   Int_t i=(Int_t) ( (theta-rad40)/rad100*imax);
 541   Int_t j=(Int_t) ( phi/rad360*jmax );
 542   //Int_t i=(Int_t)( ((theta-((TMath::Pi())/4.))/((TMath::Pi())/2.))*imax );
 543   //Int_t j=(Int_t)( (phi/((TMath::Pi())*2.))*jmax );
 544   if(i<0) i=0;
 545   if(i>=imax) i=imax-1;
 546   if(j<0) j=0;
 547   if(j>=jmax) j=jmax-1;
 548
 549   s(0) = 0.0028/TMath::Sin(theta)+( rl->GetRadMatrix1() )(i,j);   // 0.0028 takes into account the beam pipe
 550   s(1) = ( rl->GetRadMatrix2() )(i,j);
 551   s(2) = ( rl->GetRadMatrix3() )(i,j);
 552   s(3) = ( rl->GetRadMatrix4() )(i,j);
 553   s(4) = ( rl->GetRadMatrix5() )(i,j);
 554   s(5) = ( rl->GetRadMatrix6() )(i,j);
 555
 556
 557   */
 558
 559   //for(int k=0; k<6; k++) cout<<s(k)<<" "; getchar();
 560
 561   //if(phi>60) {cout<<" phi = "<<phi<<"\n"; getchar();}
 562   //if(theta<45 || theta>135) {cout<<" theta = "<<theta<<"\n"; getchar();}
 563   //cout<<" dentro AddEl: phi, theta = "<<phi<<" "<<theta<<"\n"; getchar();
 564
 565   Double_t cl=1.+fX3*fX3;  // cl=1/(cosl)**2 = 1 + (tgl)**2
 566   Double_t sqcl=TMath::Sqrt(cl);
 567   Double_t pt=GetPt();
 568
 569   Double_t p2=pt*pt*cl;
 570   Double_t e=TMath::Sqrt(p2+mass*mass);
 571   Double_t beta2=p2/(p2+mass*mass);
 572
 573   Double_t dE;
 574   if(flagtot) {
 575     Double_t stot=s(0)+s(1)+s(2)+s(3)+s(4)+s(5);
 576     dE=0.153/beta2*(log(5940*beta2/(1-beta2)) - beta2)*stot*21.82*sqcl;
 577   } else {
 578     dE=0.153/beta2*(log(5940*beta2/(1-beta2)) - beta2)*s(fLayer-1)*21.82*sqcl;
 579   }
 580   dE=signdE*dE/1000.;
 581
 582   e+=dE;
 583   Double_t p=TMath::Sqrt(e*e-mass*mass);
 584   Double_t sign=1.;
 585   if(fX4 < 0.) sign=-1.;
 586   pt=sign*p/sqcl;
 587   Double_t cc=(0.299792458*0.2*fFieldFactor)/(pt*100.);
 588   fX4=cc;
 589
 590 }
 591
 592 void  AliITSTrackV1::Correct(Double_t rk) {
 593 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 594 // correct track to take into account real geometry detector
 595
 596   Double_t duepi=2.*TMath::Pi();
 597   Double_t rkm1=frtrack;
 598   Double_t aAk=ArgA(rk), aAkm1=ArgA(rkm1);
 599   Double_t ak=Arga(rk), akm1=Arga(rkm1);
 600
 601   fX0+=TMath::ASin(aAk)-TMath::ASin(aAkm1);
 602   if(fX0>duepi) fX0-=duepi;
 603   if(fX0<0.) fX0+=duepi;
 604
 605   Double_t tgl=fX3;
 606   Double_t c=fX4;
 607   Double_t cy=c/2;
 608   fX1+=tgl/cy*(TMath::ASin(ak)-TMath::ASin(akm1));
 609   frtrack=rk;
 610
 611 }
 612
 613 void AliITSTrackV1::AddMS(AliITSRad *rl, Double_t mass) {
 614 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 615 //////////   Modification of the covariance matrix to take into account multiple scattering  ///////////
 616
 617   TVector s(6);
 618
 619   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;
 620  //  s(0)=0.0026+0.00283*2.; s(1)=0.018*2.; s(2)=0.0094; s(3)=0.0095; s(4)=0.0091; s(5)=0.0087;
 621 //0.00277 is added in the first layer to take into account the energy loss in the beam pipe
 622
 623   Int_t k;
 624   //for(k=0; k<6; k++) s(k)=s(k)*1.6;   // forint
 625   for(k=0; k<6; k++) s(k)=s(k)*1.7;   // forint
 626
 627   Double_t phi=fX0;
 628  if(phi<0.174 ) s(5)=s(5)+0.012; //aggiunta provvisoria
 629   if(phi>6.1 ) s(5)=s(5)+0.012; //aggiunta provvisoria
 630    if(phi>2.96 && phi< 3.31) s(5)=s(5)+0.012; //aggiunta provvisoria
 631
 632
 633   Double_t tgl=fX3;
 634   /*
 635   Double_t theta=((TMath::Pi())/2.)-TMath::ATan(tgl);
 636   Double_t rad40=(TMath::Pi())*40./180.;      // rivedere
 637   Double_t rad100=(TMath::Pi())*100/180;
 638   Double_t rad360=(TMath::Pi())*2.;
 639   Int_t imax=rl->Getimax();
 640   Int_t jmax=rl->Getjmax();
 641   Int_t i=(Int_t) ( (theta-rad40)/rad100*imax);
 642   Int_t j=(Int_t) ( phi/rad360*jmax);
 643
 644   if(i<0) i=0;
 645   if(i>=imax) i=imax-1;
 646   if(j<0) j=0;
 647   if(j>=jmax) j=jmax-1;
 648
 649   s(0) = 0.0028/TMath::Sin(theta)+( rl->GetRadMatrix1() )(i,j);   // 0.0028 takes into account the beam pipe
 650   s(1) = ( rl->GetRadMatrix2() )(i,j);
 651   s(2) = ( rl->GetRadMatrix3() )(i,j);
 652   s(3) = ( rl->GetRadMatrix4() )(i,j);
 653   s(4) = ( rl->GetRadMatrix5() )(i,j);
 654   s(5) = ( rl->GetRadMatrix6() )(i,j);
 655    */
 656   //Double_t mass=0.1396;
 657   Int_t layer=(Int_t)GetLayer();
 658
 659   Double_t cosl=TMath::Cos(TMath::ATan(tgl));
 660   Double_t d=fX2;
 661   Double_t c=fX4;
 662   Double_t cy=c/2.;
 663   Double_t q20=1./(cosl*cosl);
 664   Double_t q30=c*tgl;
 665
 666   Double_t q40=cy*(frtrack*frtrack-d*d)/(1.+ 2.*cy*d);
 667   Double_t dd=d+cy*d*d-cy*frtrack*frtrack;
 668   Double_t dprova=frtrack*frtrack - dd*dd;
 669   Double_t q41=0.;
 670   if(dprova>0.) q41=-1./cosl*TMath::Sqrt(dprova)/(1.+ 2.*cy*d);
 671
 672   Double_t p2=(GetPt()*GetPt())/(cosl*cosl);
 673   Double_t beta2=p2/(p2+mass*mass);
 674   Double_t theta2=14.1*14.1/(beta2*p2*1.e6)*(s(layer-1)/cosl);
 675
 676   fC22+=theta2*(q40*q40+q41*q41);
 677   fC32+=theta2*q20*q40;
 678   fC33+=theta2*q20*q20;
 679   fC42+=theta2*q30*q40;
 680   fC43+=theta2*q30*q20;
 681   fC44+=theta2*q30*q30;
 682
 683 }
 684 void AliITSTrackV1::PrimaryTrack(AliITSRad *rl) {
 685 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 686 // calculation of part of covariance matrix for vertex constraint
 687
 688   Double_t rlayer[6];
 689
 690   rlayer[0]=4.; rlayer[1]=7.;  rlayer[2]=14.9;  rlayer[3]=23.8;
 691   rlayer[4]=39.1;  rlayer[5]=43.6;
 692
 693   Double_t cy=fX4/2.;
 694   Double_t tgl=(fX1-fZv)*cy/TMath::ASin(cy*frtrack);
 695   Double_t frtrack=1.;
 696   fX0=0.;
 697   fX1=frtrack*tgl;
 698   fX2=fDv;
 699   fX3=tgl;
 700
 701   fC00=fC10=fC11=fC20=fC21=fC22=fC30=fC31=fC32=fC33=fC40=fC41=fC42=fC43=0.;
 702
 703   AddEL(rl,1.,1);
 704   fLayer=0;
 705   Int_t i;
 706   for (i=0; i<6; i++) {
 707     Propagation(rlayer[i]);
 708     fLayer++;
 709     fd2(i)=fC22;
 710     ftgl2(i)=fC33;
 711     fdtgl(i)=fC32;
 712     AddMS(rl);
 713     AddEL(rl,-1,0);
 714   }
 715 }
 716
 717
 718 Int_t AliITSTrackV1::DoNotCross(Double_t rk) const{
 719 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 720 // determine if the track cross a layer
 721
 722   Double_t c=fX4;
 723   Double_t d=fX2;
 724   Double_t cy=c/2.;
 725   return (TMath::Abs((cy*rk+(1.+cy*d)*d/rk)/(1.+2.*cy*d))>=1.)?1:0;
 726 }
 727
 728
 729 Double_t AliITSTrackV1::ArgA(Double_t rk) const {
 730 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 731 // quantity usefull in Propagation
 732   Double_t c=fX4;
 733   Double_t d=fX2;
 734   Double_t cy=c/2.;
 735   Double_t arg=(cy*rk + (1 + cy*d)*d/rk)/(1.+ 2.*cy*d);
 736   if (TMath::Abs(arg) < 1.) return arg;
 737   //cout<<"class AliITSTrack: ArgA out of range !\n";/* getchar();*/
 738   return (arg>0) ? 0.99999999999 : -0.9999999999;
 739 }
 740
 741 Double_t AliITSTrackV1::Arga(Double_t rk) const {
 742 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 743 // quantity usefull in Propagation
 744   Double_t c=fX4;
 745   Double_t d=fX2;
 746   Double_t cy=c/2.;
 747   Double_t arg=(rk*rk - d*d)/(1.+ 2.*cy*d);
 748   if (arg<0.) {/*cout<<"class AliITSTrack: Arga out of range !\n";*/ arg=0.;}
 749   return cy*TMath::Sqrt(arg);
 750 }
 751
 752 Double_t AliITSTrackV1::ArgB(Double_t rk) const {
 753 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 754 // quantity usefull in Propagation
 755   Double_t c=fX4;
 756   Double_t d=fX2;
 757   Double_t cy=c/2.;
 758   return (rk*rk - d*d)/(rk*(1.+ 2.*cy*d)*(1.+ 2.*cy*d));
 759 }
 760
 761 Double_t AliITSTrackV1::ArgC(Double_t rk) const {
 762 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 763 // quantity usefull in Propagation
 764   Double_t c=fX4;
 765   Double_t d=fX2;
 766   Double_t cy=c/2.;
 767   return  (1./rk - 2.*cy*ArgA(rk)/(1.+ 2.*cy*d));
 768 }
 769
 770
 771 Double_t AliITSTrackV1::GetPredChi2(Double_t m[2], Double_t sigma[2] ) const {
 772 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 773 // This function calculates a predicted chi2 increment.
 774 /*
 775   Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
 776   r00+=fC00; r01+=fC10; r11+=fC11;
 777
 778   Double_t det=r00*r11 - r01*r01;
 779   if (TMath::Abs(det) < 1.e-10) {
 780     Int_t n=GetNumberOfClusters();
 781     if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
 782     return 1e10;
 783   }
 784   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
 785
 786   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
 787
 788   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
 789   */
 790
 791
 792   Double_t r00=sigma[0], r01=0., r11=sigma[1];
 793   r00+=fC00; r01+=fC10; r11+=fC11;
 794   Double_t det=r00*r11-r01*r01;
 795   if(TMath::Abs(det) < 1.e-15) {cout<<" Problems on matrix in GetPredChi2 "<<det<<"\n";
 796   return 1e10;}
 797   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
 798   Double_t dphi=m[0]-fX0;
 799   Double_t dz=m[1]-fX1;
 800   Double_t chi2 = (dphi*r00*dphi +2.*r01*dphi*dz + dz*r11*dz)/det;
 801   return chi2;
 802
 803
 804
 805
 806 }