TPC/AliTPCtrackerParam.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 */
  19
  20
  21  /**************************************************************************
  22  *                                                                        *
  23  * This class builds AliTPCtrack objects from generated tracks to feed    *
  24  * ITS tracking (V2). The AliTPCtrack is built from its first hit in      *
  25  * the TPC. The track is assigned a Kalman-like covariance matrix         *
  26  * depending on its pT and pseudorapidity and track parameters are        *
  27  * smeared according to this covariance matrix.                           *
  28  * Output file contains sorted tracks, ready for matching with ITS        *
  29  *                                                                        *
  30  * For details:                                                           *
  31  * http://www.pd.infn.it/alipd/talks/soft/adIII02/TPCtrackingParam.htm    *
  32  *                                                                        *
  33  * Test macro is: AliBarrelRec_TPCparam.C                                 *
  34  *                                                                        *
  35  *  Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it    *
  36  *                                                                        *
  37  **************************************************************************/
  38 #include "AliTPCtrackerParam.h"
  39 #include "alles.h"
  40 #include "AliMagF.h"
  41 #include "AliTPCtrack.h"
  42 #include "TMatrixD.h"
  43 #include "AliKalmanTrack.h"
  44 #include "AliMagFCM.h"
  45 #include "AliGausCorr.h"
  46
  47
  48 ClassImp(AliTPCtrackerParam)
  49
  50 //-----------------------------------------------------------------
  51 AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz)
  52 {
  53 //-----------------------------------------------------------------
  54 // This is the class conctructor
  55 //-----------------------------------------------------------------
  56
  57   fColl = coll; // collision code (0: PbPb6000)
  58   fBz = Bz;     // value of the z component of L3 field (Tesla)
  59
  60 }
  61 //-----------------------------------------------------------------
  62 AliTPCtrackerParam::~AliTPCtrackerParam()
  63 {}
  64 //-----------------------------------------------------------------
  65
  66 Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out, Int_t n)
  67 {
  68 //-----------------------------------------------------------------
  69 // This function creates the TPC parameterized tracks
  70 //-----------------------------------------------------------------
  71
  72   if(fColl!=0) {
  73     cerr<<"AliTPCtrackerParam::BuildTPCtracks:  Invalid collision!\n";
  74     cerr<<"      Available:  0   ->   PbPb6000"<<endl; return 0;
  75   }
  76   if(fBz!=0.4) {
  77     cerr<<"AliTPCtrackerParam::BuildTPCtracks:  Invalid field!\n";
  78     cerr<<"      Available:  0.4"<<endl; return 0;
  79   }
  80
  81   TFile *infile=(TFile*)inp;
  82
  83   // Get gAlice object from file
  84   if(!(gAlice=(AliRun*)infile->Get("gAlice"))) {
  85     cerr<<"gAlice has not been found on galice.root !\n";
  86     return 1;
  87   }
  88
  89   AliMagFCM *fiel = (AliMagFCM*)gAlice->Field();
  90   Double_t fieval=(Double_t)fiel->SolenoidField()/10.;
  91   printf("Magnetic field is %6.2f Tesla\n",fieval);
  92   if(fBz!=fieval) {
  93     cerr<<"AliTPCtrackerParam::BuildTPCtracks:  Invalid field!"<<endl;
  94     cerr<<"Field selected is: "<<fBz<<" T\n";
  95     cerr<<"Field found on file is: "<<fieval<<" T\n";
  96     return 0;
  97   }
  98
  99   AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
 100
 101
 102   // loop over first n events in file
 103   for(Int_t evt=0; evt<n; evt++){
 104     cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
 105
 106     AliTPCtrack *tpctrack=0;
 107
 108     // tree for TPC tracks
 109     Char_t tname[100];
 110     sprintf(tname,"TreeT_TPC_%d",evt);
 111     TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
 112     tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);
 113
 114     // array for TPC tracks
 115     TObjArray tarray(20000);
 116
 117
 118     Int_t nparticles=gAlice->GetEvent(evt);
 119
 120     Bool_t done[500000];
 121     for(Int_t l=0; l<500000; l++) { done[l]=kFALSE; }
 122
 123
 124     // Get TPC detector
 125     AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
 126     Int_t ver = TPC->IsVersion();
 127     cerr<<"+++ TPC version "<<ver<<" has been found !\n";
 128     AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
 129     if(digp){
 130       delete digp;
 131       digp = new AliTPCParamSR();
 132     }
 133     else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60");
 134
 135     if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
 136     TPC->SetParam(digp);
 137
 138     // Get TreeH with hits
 139     TTree *TH=gAlice->TreeH();
 140     Int_t ntracks=(Int_t)TH->GetEntries();
 141     cerr<<"+++\n+++ Number of particles in event "<<evt<<":  "<<nparticles<<"\n+++\n+++ Number of \"primary tracks\" (entries in TreeH): "<<ntracks<<"\n+++\n\n";
 142
 143     TParticle* Part;
 144     AliTPChit* tpcHit;
 145     Double_t hPx,hPy,hPz,hPt,xg,yg,zg,xl,yl,zl;
 146     Double_t alpha;
 147     Float_t cosAlpha,sinAlpha;
 148     Int_t label,pdg,charge,bin;
 149     Int_t tracks=0;
 150     //Int_t nSel=0,nAcc=0;
 151
 152     // loop over entries in TreeH
 153     for(Int_t i=0; i<ntracks; i++) {
 154       if(i%1000==0) cerr<<"  --- Processing primary track "<<i<<" of "<<ntracks<<" ---\r";
 155       TPC->ResetHits();
 156       TH->GetEvent(i);
 157       // Get FirstHit
 158       tpcHit=(AliTPChit*)TPC->FirstHit(-1);
 159       for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) {
 160         if(tpcHit->fQ !=0.) continue;
 161         // Get particle momentum at hit
 162         hPx=tpcHit->X(); hPy=tpcHit->Y(); hPz=tpcHit->Z();
 163         hPt=TMath::Sqrt(hPx*hPx+hPy*hPy);
 164         // reject hits with Pt<mag*0.45 GeV/c
 165         if(hPt<(fBz*0.45)) continue;
 166
 167         // Get track label
 168         label=tpcHit->Track();
 169         // check if this track has already been processed
 170         if(done[label]) continue;
 171         // electric charge
 172         Part = gAlice->Particle(label);
 173         pdg = Part->GetPdgCode();
 174         if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
 175         else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
 176         else continue;
 177
 178
 179         if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break;
 180         if(tpcHit->fQ != 0.) continue;
 181         // Get global coordinates of hit
 182         xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
 183         if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;
 184
 185         // Get TPC sector, Alpha angle and local coordinates
 186         // printf("Sector %d\n",tpcHit->fSector);
 187         digp->AdjustCosSin(tpcHit->fSector,cosAlpha,sinAlpha);
 188         alpha = TMath::ATan2(sinAlpha,cosAlpha);
 189         xl = xg*cosAlpha + yg*sinAlpha;
 190         yl =-xg*sinAlpha + yg*cosAlpha;
 191         zl = zg;
 192         //printf("Alpha %f   xl %f  yl %f  zl %f\n",Alpha,xl,yl,zl);
 193
 194         // reject tracks which are not in the TPC acceptance
 195         if(TMath::Abs(zl+(244.-xl)*hPz/hPt)>252.) continue;
 196
 197         // Get bin in pT,eta
 198         bin = GetBin(hPt,Part->Eta());
 199
 200         // Apply selection according to TPC efficiency
 201         //if(TMath::Abs(pdg)==211) nAcc++;
 202         if(!SelectedTrack(pdg,hPt,Part->Eta())) continue;
 203         //if(TMath::Abs(pdg)==211) nSel++;
 204
 205         // Mark track as "done"
 206         done[label]=kTRUE; tracks++;
 207
 208         // create AliTPCtrack object
 209         tpctrack = BuildTrack(alpha,xl,yl,zl,hPx,hPy,hPz,hPt,charge,label);
 210
 211         // put track in array
 212         tarray.AddLast(tpctrack);
 213
 214       }
 215
 216     } // loop over entries in TreeH
 217
 218     TObjArray newtarray(20000);
 219
 220     // assing covariance matrixes and smear track parameters
 221     CookTracks(tarray,newtarray);
 222
 223     // sort array with TPC tracks (decreasing pT)
 224     newtarray.Sort();
 225
 226
 227     Int_t arrentr = newtarray.GetEntriesFast();
 228     //printf("\n  %d  \n\n",arrentr);
 229     for(Int_t l=0; l<arrentr; l++) {
 230       tpctrack=(AliTPCtrack*)newtarray.UncheckedAt(l);
 231       tracktree->Fill();
 232     }
 233
 234     // write the tree with tracks in the output file
 235     out->cd();
 236     tracktree->Write();
 237
 238     delete tracktree;
 239
 240     printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
 241     //printf("Average Eff: %f\n",(Float_t)nSel/nAcc);
 242
 243   } // loop on events
 244
 245   return 0;
 246 }
 247
 248 //-----------------------------------------------------------------
 249 AliTPCtrack* AliTPCtrackerParam::BuildTrack(Double_t alpha,Double_t x,
 250                                             Double_t y,Double_t z,Double_t px,
 251                                             Double_t py,Double_t pz,Double_t pt,
 252                                             Int_t ch,Int_t lab) const
 253 {
 254 //-----------------------------------------------------------------
 255 // This function uses GEANT info to set true track parameters
 256 //-----------------------------------------------------------------
 257   Double_t xref = x;
 258   Double_t xx[5],cc[15];
 259   cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.;
 260   cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
 261
 262   // Magnetic field
 263   TVector3 bfield(0.,0.,fBz);
 264
 265
 266   // radius [cm] of track projection in (x,y)
 267   Double_t rho = pt*100./0.299792458/bfield.Z();
 268   // center of track projection in local reference frame
 269   TVector3 hmom,hpos;
 270
 271
 272   // position (local) and momentum (local) at the hit
 273   // in the bending plane (z=0)
 274   hpos.SetXYZ(x,y,0.);
 275   hmom.SetXYZ(px*TMath::Cos(alpha)+py*TMath::Sin(alpha),-px*TMath::Sin(alpha)+py*TMath::Cos(alpha),0.);
 276   TVector3 vrho = hmom.Cross(bfield);
 277   vrho *= ch;
 278   vrho.SetMag(rho);
 279
 280   TVector3 vcenter = hpos+vrho;
 281
 282   Double_t x0 = vcenter.X();
 283
 284   // fX     = xref         X-coordinate of this track (reference plane)
 285   // fAlpha = Alpha        Rotation angle the local (TPC sector)
 286   // fP0    = YL           Y-coordinate of a track
 287   // fP1    = ZG           Z-coordinate of a track
 288   // fP2    = C*x0         x0 is center x in rotated frame
 289   // fP3    = Tgl          tangent of the track momentum dip angle
 290   // fP4    = C            track curvature
 291   xx[0] = y;
 292   xx[1] = z;
 293   xx[3] = pz/pt;
 294   xx[4] = -ch/rho;
 295   xx[2] = xx[4]*x0;
 296
 297   // create the object AliTPCtrack
 298   AliTPCtrack* track = new AliTPCtrack(0,xx,cc,xref,alpha);
 299   // set the label
 300   track->SetLabel(lab);
 301
 302   return track;
 303 }
 304
 305 //-----------------------------------------------------------------
 306 Bool_t AliTPCtrackerParam::SelectedTrack(Int_t pdg,Double_t pt,Double_t eta)
 307  const {
 308 //-----------------------------------------------------------------
 309 // This function makes a selection according to TPC tracking efficiency
 310 //-----------------------------------------------------------------
 311
 312   Double_t eff=0.;
 313
 314   //eff computed with | zl+(244-xl)*pz/pt | < 252
 315   Double_t effPi[27] = {0.724587,0.743389,0.619273,0.798477,0.812036,0.823195,0.771437,0.775826,0.784136,0.809071,0.762001,0.774576,0.848834,0.787201,0.792548,0.942089,0.951631,0.951085,0.960885,0.971451,0.969103,0.983245,0.978939,0.988706,0.990852,0.985679,0.993606};
 316   Double_t effK[18]  = {0.377934,0.363962,0.321721,0.518784,0.547459,0.517878,0.612704,0.619101,0.620894,0.733411,0.732128,0.750373,0.790630,0.806565,0.791353,0.967486,0.970483,0.974527};
 317   Double_t effP[15]  = {0.131173,0.165114,0.229658,0.365357,0.412989,0.483297,0.454614,0.505173,0.658615,0.694753,0.730661,0.815680,0.873461,0.887227,0.899324};
 318   Double_t effEl[15] = {0.835549,0.853746,0.718207,0.835230,0.831489,0.862222,0.757783,0.747301,0.824096,0.867949,0.871891,0.808480,0.890625,0.911765,0.973684};
 319   Double_t effMu[15] = {0.553486,0.641392,0.609932,0.591126,0.706729,0.750755,0.747952,0.729051,0.760849,0.898810,0.737500,0.830357,0.735294,0.800000,0.882353};
 320
 321
 322   if(TMath::Abs(pdg)==211)  eff = LinearInterpolation(9,effPi,pt,eta);
 323   if(TMath::Abs(pdg)==321)  eff = LinearInterpolation(6,effK,pt,eta);
 324   if(TMath::Abs(pdg)==2212) eff = LinearInterpolation(5,effP,pt,eta);
 325   if(TMath::Abs(pdg)==11)   eff = LinearInterpolation(5,effEl,pt,eta);
 326   if(TMath::Abs(pdg)==13)   eff = LinearInterpolation(5,effMu,pt,eta);
 327
 328   if(gRandom->Rndm() < eff) return kTRUE;
 329
 330   return kFALSE;
 331 }
 332
 333 //-----------------------------------------------------------------
 334 Double_t AliTPCtrackerParam::LinearInterpolation(Int_t ptBins,Double_t *value,
 335                                                  Double_t trkPt,Double_t trkEta) const
 336 {
 337 //-----------------------------------------------------------------
 338 // This function makes a linear interpolation
 339 //-----------------------------------------------------------------
 340   Double_t intValue=0,intValue1=0,intValue2=0;
 341   Int_t etaSide = (TMath::Abs(trkEta)<.45 ? 0 : 1);
 342   Double_t eta[3]={0.15,0.45,0.75};
 343   Double_t pt[9]={0.244,0.390,0.676,1.190,2.36,4.,6.,10.,20.};
 344   if(ptBins==6) pt[5]=10.;
 345
 346   for(Int_t i=0; i<ptBins; i++) {
 347     if(trkPt<pt[i]) {
 348       if(i==0) i=1;
 349       intValue1 = value[3*i-3+etaSide]+(value[3*i-3+etaSide]-value[3*i+etaSide])/(pt[i-1]-pt[i])*(trkPt-pt[i-1]);
 350       intValue2 = value[3*i-3+etaSide+1]+(value[3*i-3+etaSide+1]-value[3*i+etaSide+1])/(pt[i-1]-pt[i])*(trkPt-pt[i-1]);
 351
 352       intValue = intValue1+(intValue1-intValue2)/(eta[etaSide]-eta[etaSide+1])*(trkEta-eta[etaSide]);
 353       break;
 354     }
 355     if(i==ptBins-1) {
 356       i=ptBins-2;
 357       intValue1 = value[3*i-3+etaSide]+(value[3*i-3+etaSide]-value[3*i+etaSide])/(pt[i-1]-pt[i])*(trkPt-pt[i-1]);
 358       intValue2 = value[3*i-3+etaSide+1]+(value[3*i-3+etaSide+1]-value[3*i+etaSide+1])/(pt[i-1]-pt[i])*(trkPt-pt[i-1]);
 359
 360       intValue = intValue1+(intValue1-intValue2)/(eta[etaSide]-eta[etaSide+1])*(trkEta-eta[etaSide]);
 361       break;
 362     }
 363   }
 364
 365   return intValue;
 366 }
 367
 368 //-----------------------------------------------------------------
 369 Int_t AliTPCtrackerParam::GetBin(Double_t pt,Double_t eta) const {
 370 //-----------------------------------------------------------------
 371 // This function tells bin number in a grid (pT,eta)
 372 //-----------------------------------------------------------------
 373   if(TMath::Abs(eta)<0.3) {
 374     if(pt<0.3)            return 0;
 375     if(pt>=0.3 && pt<0.5) return 3;
 376     if(pt>=0.5 && pt<1.)  return 6;
 377     if(pt>=1. && pt<1.5)  return 9;
 378     if(pt>=1.5 && pt<3.)  return 12;
 379     if(pt>=3. && pt<5.)   return 15;
 380     if(pt>=5. && pt<7.)   return 18;
 381     if(pt>=7. && pt<15.)  return 21;
 382     if(pt>=15.)           return 24;
 383   }
 384   if(TMath::Abs(eta)>=0.3 && TMath::Abs(eta)<0.6) {
 385     if(pt<0.3)            return 1;
 386     if(pt>=0.3 && pt<0.5) return 4;
 387     if(pt>=0.5 && pt<1.)  return 7;
 388     if(pt>=1. && pt<1.5)  return 10;
 389     if(pt>=1.5 && pt<3.)  return 13;
 390     if(pt>=3. && pt<5.)   return 16;
 391     if(pt>=5. && pt<7.)   return 19;
 392     if(pt>=7. && pt<15.)  return 22;
 393     if(pt>=15.)           return 25;
 394   }
 395   if(TMath::Abs(eta)>=0.6) {
 396     if(pt<0.3)            return 2;
 397     if(pt>=0.3 && pt<0.5) return 5;
 398     if(pt>=0.5 && pt<1.)  return 8;
 399     if(pt>=1. && pt<1.5)  return 11;
 400     if(pt>=1.5 && pt<3.)  return 14;
 401     if(pt>=3. && pt<5.)   return 17;
 402     if(pt>=5. && pt<7.)   return 20;
 403     if(pt>=7. && pt<15.)  return 23;
 404     if(pt>=15.)           return 26;
 405   }
 406
 407   return -1;
 408
 409 }
 410
 411 //-----------------------------------------------------------------
 412 TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,Double_t eta) const
 413 {
 414 //-----------------------------------------------------------------
 415 // This function stretches the covariance matrix according to the pulls
 416 //-----------------------------------------------------------------
 417   TMatrixD covMat(5,5);
 418
 419   covMat(0,0)=cc[0];
 420   covMat(1,0)=cc[1];  covMat(0,1)=covMat(1,0);
 421   covMat(1,1)=cc[2];
 422   covMat(2,0)=cc[3];  covMat(0,2)=covMat(2,0);
 423   covMat(2,1)=cc[4];  covMat(1,2)=covMat(2,1);
 424   covMat(2,2)=cc[5];
 425   covMat(3,0)=cc[6];  covMat(0,3)=covMat(3,0);
 426   covMat(3,1)=cc[7];  covMat(1,3)=covMat(3,1);
 427   covMat(3,2)=cc[8];  covMat(2,3)=covMat(3,2);
 428   covMat(3,3)=cc[9];
 429   covMat(4,0)=cc[10]; covMat(0,4)=covMat(4,0);
 430   covMat(4,1)=cc[11]; covMat(1,4)=covMat(4,1);
 431   covMat(4,2)=cc[12]; covMat(2,4)=covMat(4,2);
 432   covMat(4,3)=cc[13]; covMat(3,4)=covMat(4,3);
 433   covMat(4,4)=cc[14];
 434
 435   TMatrixD stretchMat(5,5);
 436   for(Int_t k=0;k<5;k++) {
 437     for(Int_t l=0;l<5;l++) {
 438       stretchMat(k,l)=0.;
 439     }
 440   }
 441
 442
 443   Double_t s00[27]={1.3553,1.2973,1.2446,1.3428,1.3031,1.2345,1.3244,1.2681,1.2046,1.3046,1.2430,1.1606,1.2462,1.2104,1.1082,1.2207,1.1189,1.0789,1.2079,1.1300,1.0426,1.1502,1.1122,1.0559,1.1419,1.1072,1.0208};
 444   Double_t s11[27]={1.0890,1.1463,1.2313,1.1149,1.1597,1.2280,1.1225,1.1584,1.2329,1.1149,1.1550,1.2369,1.1095,1.1353,1.2050,1.0649,1.0858,1.1546,1.0663,1.0672,1.1340,1.0416,1.0738,1.0945,1.0663,1.0654,1.0909};
 445   Double_t s22[27]={1.1709,1.1367,1.0932,1.2247,1.1832,1.1247,1.2470,1.2017,1.1441,1.2202,1.1653,1.1050,1.1819,1.1269,1.0583,1.1546,1.0621,0.9928,1.1305,1.0512,0.9576,1.0995,1.0445,0.9884,1.0968,1.0368,0.9574};
 446   Double_t s33[27]={0.9720,0.9869,1.0271,1.0030,1.0223,1.0479,1.0164,1.0305,1.0559,1.0339,1.0450,1.0686,1.0450,1.0507,1.0784,1.0334,1.0208,1.0863,1.0252,1.0114,1.0835,0.9854,1.0144,1.0507,1.0124,1.0159,1.0464};
 447   Double_t s44[27]={1.1104,1.0789,1.0479,1.1597,1.1234,1.0728,1.2087,1.1602,1.1041,1.1942,1.1428,1.0831,1.1572,1.1036,1.0431,1.1296,1.0498,0.9844,1.1145,1.0266,0.9489,1.0959,1.0450,0.9875,1.0775,1.0266,0.9406};
 448
 449
 450   stretchMat(0,0) = LinearInterpolation(9,s00,pt,eta);
 451   stretchMat(1,1) = LinearInterpolation(9,s11,pt,eta);
 452   stretchMat(2,2) = LinearInterpolation(9,s22,pt,eta);
 453   stretchMat(3,3) = LinearInterpolation(9,s33,pt,eta);
 454   stretchMat(4,4) = LinearInterpolation(9,s44,pt,eta);
 455
 456   TMatrixD mat(stretchMat,TMatrixD::kMult,covMat);
 457   TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat);
 458
 459
 460   return covMatSmear;
 461 }
 462
 463 //-----------------------------------------------------------------
 464 void AliTPCtrackerParam::SmearTrack(Double_t* xx,Double_t* xxsm,TMatrixD cov)  const {
 465 //-----------------------------------------------------------------
 466 // This function smears track parameters according to streched cov. matrix
 467 //-----------------------------------------------------------------
 468   AliGausCorr *corgen = new AliGausCorr(cov,5);
 469   TArrayD corr(5);
 470   corgen->GetGaussN(corr);
 471   delete corgen;
 472   corgen = 0;
 473
 474   for(Int_t l=0;l<5;l++) {
 475     xxsm[l] = xx[l]+corr[l];
 476   }
 477
 478   return;
 479 }
 480
 481 //-----------------------------------------------------------------
 482 void AliTPCtrackerParam::CookTracks(TObjArray& tarray,TObjArray& newtarray)
 483 const {
 484 //-----------------------------------------------------------------
 485 // This function deals with covariance matrix and smearing
 486 //-----------------------------------------------------------------
 487
 488   TString s("$ALICE_ROOT/TPC/CovMatrixDB_");
 489   if(fColl==0) s.Append("PbPb6000");
 490   if(fBz==0.4) s.Append("_B0.4T.root");
 491   // open file with matrixes DB
 492   TFile* fileDB = TFile::Open(s.Data());
 493
 494   AliTPCtrack* track = 0;
 495   Int_t arrayEntries = (Int_t)tarray.GetEntriesFast();
 496
 497   for(Int_t k=0; k<arrayEntries; k++) {
 498     track=(AliTPCtrack*)tarray.UncheckedAt(k);
 499
 500     Double_t pt = 1/TMath::Abs(track->Get1Pt());
 501     Double_t eta = -TMath::Log(TMath::Tan(0.25*TMath::Pi()-0.5*TMath::ATan(track->GetTgl())));
 502
 503     Int_t bin = GetBin(pt,eta);
 504
 505     // string with tree name
 506     TString str("Tree_bin");
 507     str+=bin;
 508
 509     // get the right tree
 510     TTree* covTree = (TTree*)fileDB->Get(str.Data());
 511     TArrayF* matrix = new TArrayF;
 512     covTree->SetBranchAddress("matrixes",&matrix);
 513
 514     // get random entry from the tree
 515     Int_t treeEntries = (Int_t)covTree->GetEntries();
 516     covTree->GetEvent(gRandom->Integer(treeEntries));
 517
 518     Double_t xref,alpha,xx[5],xxsm[5],cc[15];
 519     Int_t lab;
 520     // get P and Cosl from track
 521     Double_t cosl = TMath::Cos(TMath::ATan(track->GetTgl()));
 522     Double_t p = 1/TMath::Abs(track->Get1Pt())/cosl;
 523
 524     // get covariance matrix from regularized matrix
 525     cc[0]=matrix->At(0)*(1.588e-3+1.911e-4/TMath::Power(p,1.5));
 526     cc[1]=matrix->At(1);
 527     cc[2]=matrix->At(2)*(1.195e-3+0.8102e-3/p);
 528     cc[3]=matrix->At(3)*(7.275e-5+1.181e-5/TMath::Power(p,1.5));
 529     cc[4]=matrix->At(4);
 530     cc[5]=matrix->At(5)*(5.163e-6+1.138e-6/TMath::Power(p,2.)/cosl);
 531     cc[6]=matrix->At(6);
 532     cc[7]=matrix->At(7)*(1.176e-5+1.175e-5/TMath::Power(p,1.5)/cosl/cosl/cosl);
 533     cc[8]=matrix->At(8);
 534     cc[9]=matrix->At(9)*(1.289e-7+4.636e-7/TMath::Power(p,1.7)/cosl/cosl/cosl/cosl);
 535     cc[10]=matrix->At(10)*(4.056e-7+4.379e-8/TMath::Power(p,1.5));
 536     cc[11]=matrix->At(11);
 537     cc[12]=matrix->At(12)*(3.049e-8+8.244e-9/TMath::Power(p,2.)/cosl);
 538     cc[13]=matrix->At(13);
 539     cc[14]=matrix->At(14)*(1.847e-10+5.822e-11/TMath::Power(p,2.)/cosl/cosl/cosl);
 540
 541     TMatrixD covMatSmear(5,5);
 542
 543     covMatSmear = GetSmearingMatrix(cc,pt,eta);
 544
 545     // get track original parameters
 546     xref=track->GetX();
 547     alpha=track->GetAlpha();
 548     xx[0]=track->GetY();
 549     xx[1]=track->GetZ();
 550     xx[2]=track->GetX()*track->GetC()-track->GetSnp();
 551     xx[3]=track->GetTgl();
 552     xx[4]=track->GetC();
 553     lab=track->GetLabel();
 554
 555     // use smearing matrix to smear the original parameters
 556     SmearTrack(xx,xxsm,covMatSmear);
 557
 558     AliTPCtrack* tpctrack = new AliTPCtrack(0,xxsm,cc,xref,alpha);
 559     tpctrack->SetLabel(lab);
 560
 561     // fill the array
 562     newtarray.AddLast(tpctrack);
 563
 564     delete matrix;
 565
 566   }
 567
 568   fileDB->Close();
 569   delete fileDB;
 570
 571
 572   return;
 573 }
 574 //-----------------------------------------------------------------
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586