[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerParam.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
/*
$Log$
*/
/*************************************************************************
 *                                                                        *
 * This class builds AliTPCtrack objects from generated tracks to feed    *
 * ITS tracking (V2). The AliTPCtrack is built from its first hit in      *
 * the TPC. The track is assigned a Kalman-like covariance matrix         *
 * depending on its pT and pseudorapidity and track parameters are        *
 * smeared according this covariance matrix.                              *
 * Output file contains sorted tracks, ready for matching with ITS        *
 *                                                                        *
 * Test macro is: AliBarrelRec_TPCparam.C                                 * 
 *                                                                        *
 *  Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it    * 
 *                                                                        *
 **************************************************************************/
#include "AliTPCtrackerParam.h"
#include "alles.h"
#include "AliMagF.h"
#include "AliTPCtrack.h"
#include "TMatrixD.h"
#include "AliKalmanTrack.h"

ClassImp(AliTPCtrackerParam)

//-----------------------------------------------------------------
AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz) 
{
  fColl = coll; // collision code (0: PbPb6000)
  fBz = Bz;     // value of the z component of L3 field (Tesla)
  
}
//-----------------------------------------------------------------
AliTPCtrackerParam::~AliTPCtrackerParam() 
{}
//-----------------------------------------------------------------

Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out, Int_t n)
{
 
  if(fColl!=0) { 
    cerr<<"AliTPCtrackerParam::BuildTPCtracks:  Invalid collision!\n
                 Available:  0   ->   PbPb6000"<<endl; return 0; 
  }
  if(fBz!=0.4) {
    cerr<<"AliTPCtrackerParam::BuildTPCtracks:  Invalid field!\n
                 Available:  0.4"<<endl; return 0;
  }

  TFile *infile=(TFile*)inp;

  AliKalmanTrack::SetConvConst(100/0.299792458/fBz);

  // Get gAlice object from file
  if(!(gAlice=(AliRun*)infile->Get("gAlice"))) {
    cerr<<"gAlice has not been found on galice.root !\n";
    return 1;
  }

  // Set random number generator seed
  TDatime t;
  gRandom->SetSeed(t.Get());


  // loop over first n events in file
  for(Int_t evt=0; evt<n; evt++){
    cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";

    AliTPCtrack *tpctrack=0;

    // tree for TPC tracks
    Char_t tname[100];
    sprintf(tname,"TreeT_TPC_%d",evt);
    TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
    tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);

    // array for TPC tracks
    TObjArray tarray(20000);


    Int_t nparticles=gAlice->GetEvent(evt);   

    Bool_t done[500000];
    for(Int_t l=0; l<500000; l++) { done[l]=kFALSE; }


    // Get TPC detector 
    AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
    Int_t ver = TPC->IsVersion(); 
    cerr<<"+++ TPC version "<<ver<<" has been found !\n";
    AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
    if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
    TPC->SetParam(digp);

    // Get TreeH with hits
    TTree *TH=gAlice->TreeH(); 
    Int_t ntracks=(Int_t)TH->GetEntries();
    cerr<<"+++\n+++ Number of particles in event "<<evt<<":  "<<nparticles<<"\n+++\n+++ Number of \"primary tracks\" (entries in TreeH): "<<ntracks<<"\n+++\n\n";

    TParticle* Part;
    AliTPChit* tpcHit;
    Double_t hPx,hPy,hPz,hPt,xg,yg,zg,xl,yl,zl;
    Double_t Alpha;
    Float_t cosAlpha,sinAlpha;
    Int_t label,Pdg,charge,bin;
    Int_t tracks=0;
    //Int_t N1=0,N2=0;

    // loop over entries in TreeH
    for(Int_t i=0; i<ntracks; i++) {
      if(i%1000==0) cerr<<"  --- Processing primary track "<<i<<" of "<<ntracks<<" ---\r";
      TPC->ResetHits();
      TH->GetEvent(i);
      // Get FirstHit
      tpcHit=(AliTPChit*)TPC->FirstHit(-1);
      for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) {
	if(tpcHit->fQ !=0.) continue;
	// Get particle momentum at hit
	hPx=tpcHit->X(); hPy=tpcHit->Y(); hPz=tpcHit->Z();
	hPt=TMath::Sqrt(hPx*hPx+hPy*hPy);
	// reject hits with Pt<mag*0.45 GeV/c
	if(hPt<(fBz*0.45)) continue;

	// Get track label
	label=tpcHit->Track();
	// check if this track has already been processed
	if(done[label]) continue;
	// electric charge
	Part = gAlice->Particle(label);
	Pdg = Part->GetPdgCode();
	if(Pdg>200 || Pdg==-11 || Pdg==-13) { charge=1; }
	else if(Pdg<-200 || Pdg==11 || Pdg==13) { charge=-1; }
	else continue;


	if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break;
	if(tpcHit->fQ != 0.) continue;
	// Get global coordinates of hit
	xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
	if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;

	// Get TPC sector, Alpha angle and local coordinates
	// printf("Sector %d\n",tpcHit->fSector);
	digp->AdjustCosSin(tpcHit->fSector,cosAlpha,sinAlpha);
	Alpha = TMath::ATan2(sinAlpha,cosAlpha);
	xl = xg*cosAlpha + yg*sinAlpha;
	yl =-xg*sinAlpha + yg*cosAlpha;
	zl = zg;
	//printf("Alpha %f   xl %f  yl %f  zl %f\n",Alpha,xl,yl,zl);

	// reject tracks which are not in the TPC acceptance
	if(TMath::Abs(zl+(244.-xl)*hPz/hPt)>252.) continue;

	// Get bin in pT,eta
	bin = GetBin(hPt,Part->Eta());

	// Apply selection according to TPC efficiency
	//if(abs(Pdg)==211) N1++;
	if(!SelectedTrack(Pdg,hPt,Part->Eta())) continue; 
	//if(abs(Pdg)==211) N2++;

	// Mark track as "done"
	done[label]=kTRUE; tracks++;
 
	// create AliTPCtrack object
	tpctrack = BuildTrack(Alpha,xl,yl,zl,hPx,hPy,hPz,hPt,charge,label);

	// put track in array
	tarray.AddLast(tpctrack);

      }
 
    } // loop over entries in TreeH

    TObjArray newtarray(20000);

    // assing covariance matrixes and smear track parameters
    CookTracks(tarray,newtarray);

    // sort array with TPC tracks (decreasing pT)
    newtarray.Sort();


    Int_t arrentr = newtarray.GetEntriesFast();
    //printf("\n  %d  \n\n",arrentr);
    for(Int_t l=0; l<arrentr; l++) {
      tpctrack=(AliTPCtrack*)newtarray.UncheckedAt(l);
      tracktree->Fill();
    }

    // write the tree with tracks in the output file
    out->cd();
    tracktree->Write();

    delete tracktree;

    printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
    //printf("Average Eff: %f\n",(Float_t)N2/N1);

  } // loop on events

  return 0;
}

//-----------------------------------------------------------------
AliTPCtrack* AliTPCtrackerParam::BuildTrack(Double_t alpha,Double_t x,
					    Double_t y,Double_t z,Double_t Px,
					    Double_t Py,Double_t Pz,Double_t Pt,
					    Int_t ch,Int_t lab) 
{  
  Double_t xref = x;
  Double_t xx[5],cc[15];
  cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.;
  cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
  
  // Magnetic field
  TVector3 Bfield(0.,0.,fBz);
  
  
  // radius [cm] of track projection in (x,y) 
  Double_t rho = Pt*100/0.299792458/Bfield.Z();
  // center of track projection in local reference frame
  TVector3 MomH,PosH;


  // position (local) and momentum (local) at the hit
  // in the bending plane (z=0)
  PosH.SetXYZ(x,y,0.);
  MomH.SetXYZ(Px*TMath::Cos(alpha)+Py*TMath::Sin(alpha),-Px*TMath::Sin(alpha)+Py*TMath::Cos(alpha),0.);
  TVector3 Vrho = MomH.Cross(Bfield);
  Vrho *= ch;
  Vrho.SetMag(rho);

  TVector3 Vcenter = PosH+Vrho;

  Double_t x0 = Vcenter.X();

  // fX     = xref         X-coordinate of this track (reference plane)
  // fAlpha = Alpha        Rotation angle the local (TPC sector) 
  // fP0    = YL           Y-coordinate of a track
  // fP1    = ZG           Z-coordinate of a track
  // fP2    = C*x0         x0 is center x in rotated frame
  // fP3    = Tgl          tangent of the track momentum dip angle
  // fP4    = C            track curvature
  xx[0] = y;
  xx[1] = z;
  xx[3] = Pz/Pt;
  xx[4] = -ch/rho;
  xx[2] = xx[4]*x0;

  // create the object AliTPCtrack    
  AliTPCtrack* track = new AliTPCtrack(0,xx,cc,xref,alpha);
  // set the label
  track->SetLabel(lab);

  return track;
}

//-----------------------------------------------------------------
Bool_t AliTPCtrackerParam::SelectedTrack(Int_t pdg,Double_t pt,Double_t eta) 
{
  Double_t eff=0.;
  
  //eff computed with | zl+(244-xl)*pz/pt | < 252
  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};
  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};
  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};
  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};
  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};
 

  if(TMath::Abs(pdg)==211)  eff = LinearInterpolation(9,effPi,pt,eta);
  if(TMath::Abs(pdg)==321)  eff = LinearInterpolation(6,effK,pt,eta);
  if(TMath::Abs(pdg)==2212) eff = LinearInterpolation(5,effP,pt,eta);
  if(TMath::Abs(pdg)==11)   eff = LinearInterpolation(5,effEl,pt,eta);
  if(TMath::Abs(pdg)==13)   eff = LinearInterpolation(5,effMu,pt,eta);

  if(gRandom->Rndm() < eff) return kTRUE;

  return kFALSE;
}

//-----------------------------------------------------------------
Double_t AliTPCtrackerParam::LinearInterpolation(Int_t ptBins,Double_t *value,
						 Double_t pt,Double_t eta) 
{
  Double_t kValue=0,kValue1=0,kValue2=0;
  Int_t kEtaSide = (TMath::Abs(eta)<.45 ? 0 : 1);
  Double_t kEta[3]={0.15,0.45,0.75};
  Double_t kPt[9]={0.244,0.390,0.676,1.190,2.36,4.,6.,10.,20.};
  if(ptBins==6) kPt[5]=10.;

  for(Int_t i=0; i<ptBins; i++) {
    if(pt<kPt[i]) {
      if(i==0) i=1;
      kValue1 = value[3*i-3+kEtaSide]+(value[3*i-3+kEtaSide]-value[3*i+kEtaSide])/(kPt[i-1]-kPt[i])*(pt-kPt[i-1]);
      kValue2 = value[3*i-3+kEtaSide+1]+(value[3*i-3+kEtaSide+1]-value[3*i+kEtaSide+1])/(kPt[i-1]-kPt[i])*(pt-kPt[i-1]);

      kValue = kValue1+(kValue1-kValue2)/(kEta[kEtaSide]-kEta[kEtaSide+1])*(eta-kEta[kEtaSide]);
      break;
    }
    if(i==ptBins-1) { 
      i=ptBins-2;
      kValue1 = value[3*i-3+kEtaSide]+(value[3*i-3+kEtaSide]-value[3*i+kEtaSide])/(kPt[i-1]-kPt[i])*(pt-kPt[i-1]);
      kValue2 = value[3*i-3+kEtaSide+1]+(value[3*i-3+kEtaSide+1]-value[3*i+kEtaSide+1])/(kPt[i-1]-kPt[i])*(pt-kPt[i-1]);

      kValue = kValue1+(kValue1-kValue2)/(kEta[kEtaSide]-kEta[kEtaSide+1])*(eta-kEta[kEtaSide]);
      break;
    }
  }

  return kValue;
}

//-----------------------------------------------------------------
Int_t AliTPCtrackerParam::GetBin(Double_t pt,Double_t eta) {

  if(TMath::Abs(eta)<0.3) {
    if(pt<0.3)            return 0;
    if(pt>=0.3 && pt<0.5) return 3;
    if(pt>=0.5 && pt<1.)  return 6;
    if(pt>=1. && pt<1.5)  return 9;
    if(pt>=1.5 && pt<3.)  return 12;
    if(pt>=3. && pt<5.)   return 15;
    if(pt>=5. && pt<7.)   return 18;
    if(pt>=7. && pt<15.)  return 21;
    if(pt>=15.)           return 24;
  }
  if(TMath::Abs(eta)>=0.3 && TMath::Abs(eta)<0.6) {
    if(pt<0.3)            return 1;
    if(pt>=0.3 && pt<0.5) return 4;
    if(pt>=0.5 && pt<1.)  return 7;
    if(pt>=1. && pt<1.5)  return 10;
    if(pt>=1.5 && pt<3.)  return 13;
    if(pt>=3. && pt<5.)   return 16;
    if(pt>=5. && pt<7.)   return 19;
    if(pt>=7. && pt<15.)  return 22;
    if(pt>=15.)           return 25;
  }
  if(TMath::Abs(eta)>=0.6) {
    if(pt<0.3)            return 2;
    if(pt>=0.3 && pt<0.5) return 5;
    if(pt>=0.5 && pt<1.)  return 8;
    if(pt>=1. && pt<1.5)  return 11;
    if(pt>=1.5 && pt<3.)  return 14;
    if(pt>=3. && pt<5.)   return 17;
    if(pt>=5. && pt<7.)   return 20;
    if(pt>=7. && pt<15.)  return 23;
    if(pt>=15.)           return 26;
  }

  return -1;

}

//-----------------------------------------------------------------

TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,Double_t eta) 
{
  TMatrixD CovMat(5,5);

  CovMat(0,0)=cc[0];
  CovMat(1,0)=cc[1];  CovMat(0,1)=CovMat(1,0);
  CovMat(1,1)=cc[2];
  CovMat(2,0)=cc[3];  CovMat(0,2)=CovMat(2,0);
  CovMat(2,1)=cc[4];  CovMat(1,2)=CovMat(2,1);
  CovMat(2,2)=cc[5];
  CovMat(3,0)=cc[6];  CovMat(0,3)=CovMat(3,0);
  CovMat(3,1)=cc[7];  CovMat(1,3)=CovMat(3,1);
  CovMat(3,2)=cc[8];  CovMat(2,3)=CovMat(3,2);
  CovMat(3,3)=cc[9];
  CovMat(4,0)=cc[10]; CovMat(0,4)=CovMat(4,0);
  CovMat(4,1)=cc[11]; CovMat(1,4)=CovMat(4,1);
  CovMat(4,2)=cc[12]; CovMat(2,4)=CovMat(4,2);
  CovMat(4,3)=cc[13]; CovMat(3,4)=CovMat(4,3);
  CovMat(4,4)=cc[14];

  TMatrixD ScaleMat(5,5);
  for(Int_t k=0;k<5;k++) {
    for(Int_t l=0;l<5;l++) {
      ScaleMat(k,l)=0.;
    }
  }

  
  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};
  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};
  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};
  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};
  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};


  ScaleMat(0,0) = LinearInterpolation(9,s00,pt,eta); 
  ScaleMat(1,1) = LinearInterpolation(9,s11,pt,eta); 
  ScaleMat(2,2) = LinearInterpolation(9,s22,pt,eta); 
  ScaleMat(3,3) = LinearInterpolation(9,s33,pt,eta); 
  ScaleMat(4,4) = LinearInterpolation(9,s44,pt,eta); 

  TMatrixD Mat(ScaleMat,TMatrixD::kMult,CovMat);
  TMatrixD CovMatSmear(Mat,TMatrixD::kMult,ScaleMat);


  return CovMatSmear;
}

//-----------------------------------------------------------------
void AliTPCtrackerParam::SmearTrack(Double_t* xx,Double_t* xxsm,TMatrixD cov) {

  // build triangular matrix Bmat
  TMatrixD* Bmat = new TMatrixD(5,5);

  for(Int_t j=0;j<5;j++){
    Double_t accum = 0;
    for(Int_t k=0;k<j;k++){
      accum += (*Bmat)(j,k)* (*Bmat)(j,k);
    }
    (*Bmat)(j,j)=TMath::Sqrt(TMath::Abs(cov(j,j)-accum));
    for(Int_t i=j+1;i<5;i++){
      accum = 0;
      for(Int_t k=0;k<j;k++){
	accum+=(*Bmat)(i,k)* (*Bmat)(j,k);
      }
      (*Bmat)(i,j) = (cov(i,j)-accum) / (*Bmat)(j,j);
    }
  }

  // get array of numbers with normal distribution
  TArrayD norm;
  norm.Set(5);
  for(Int_t l=0;l<5;l++) {
    norm[l]=gRandom->Gaus();
  }

  // use Bmat matrix to generate correlated numbers
  TArrayD corr;
  corr.Set(5);
  for(Int_t i=0;i<5;i++){
    corr[i]=0;
    for(Int_t j=0;j<=i;j++) corr[i] += (*Bmat)(i,j)*norm[j];
  }

  delete Bmat;

  for(Int_t l=0;l<5;l++) {
    xxsm[l] = xx[l]+corr[l];
  }

  return;
}

//-----------------------------------------------------------------
void AliTPCtrackerParam::CookTracks(TObjArray& tarray,TObjArray& newtarray) 
{
  TString* s = new TString("$ALICE_ROOT/TPC/CovMatrixDB_");
  if(fColl==0) s->Append("PbPb6000");
  if(fBz==0.4) s->Append("_B0.4T.root");  
  // open file with matrixes DB
  TFile* DBfile = new TFile(s->Data());  

  AliTPCtrack* track=new AliTPCtrack;
  Int_t entr = (Int_t)tarray.GetEntriesFast();

  for(Int_t k=0; k<entr; k++) {
    track=(AliTPCtrack*)tarray.UncheckedAt(k);  

    Double_t pt = 1/TMath::Abs(track->Get1Pt());
    Double_t eta = -TMath::Log(TMath::Tan(0.25*TMath::Pi()-0.5*TMath::ATan(track->GetTgl())));

    Int_t bin = GetBin(pt,eta);

    // string with tree name
    TString str("Tree_bin");
    str+=bin;
 
    // get the right tree
    TTree* CovTree = (TTree*)DBfile->Get(str.Data());
    TArrayF* matrix = new TArrayF; 
    CovTree->SetBranchAddress("matrixes",&matrix);

    // get random entry from the tree
    Int_t Entries = (Int_t)CovTree->GetEntries();
    CovTree->GetEvent(gRandom->Integer(Entries));

    Double_t xref,alpha,xx[5],xxsm[5],cc[15];
    Int_t lab;
    // get P and Cosl from track
    Double_t Cosl = TMath::Cos(TMath::ATan(track->GetTgl()));
    Double_t P = 1/TMath::Abs(track->Get1Pt())/Cosl;
    
    // get covariance matrix from regularized matrix
    cc[0]=matrix->At(0)*(1.588e-3+1.911e-4/TMath::Power(P,1.5));
    cc[1]=matrix->At(1);
    cc[2]=matrix->At(2)*(1.195e-3+0.8102e-3/P);
    cc[3]=matrix->At(3)*(7.275e-5+1.181e-5/TMath::Power(P,1.5));
    cc[4]=matrix->At(4);
    cc[5]=matrix->At(5)*(5.163e-6+1.138e-6/TMath::Power(P,2.)/Cosl);
    cc[6]=matrix->At(6);
    cc[7]=matrix->At(7)*(1.176e-5+1.175e-5/TMath::Power(P,1.5)/Cosl/Cosl/Cosl);
    cc[8]=matrix->At(8);
    cc[9]=matrix->At(9)*(1.289e-7+4.636e-7/TMath::Power(P,1.7)/Cosl/Cosl/Cosl/Cosl);
    cc[10]=matrix->At(10)*(4.056e-7+4.379e-8/TMath::Power(P,1.5));
    cc[11]=matrix->At(11);
    cc[12]=matrix->At(12)*(3.049e-8+8.244e-9/TMath::Power(P,2.)/Cosl);
    cc[13]=matrix->At(13);
    cc[14]=matrix->At(14)*(1.847e-10+5.822e-11/TMath::Power(P,2.)/Cosl/Cosl/Cosl);
  
    TMatrixD CovMatSmear(5,5);
    
    CovMatSmear = GetSmearingMatrix(cc,pt,eta);

    // get track original parameters
    xref=track->GetX();
    alpha=track->GetAlpha();
    xx[0]=track->GetY();
    xx[1]=track->GetZ();
    xx[2]=track->GetX()*track->GetC()-track->GetSnp();
    xx[3]=track->GetTgl();
    xx[4]=track->GetC();
    lab=track->GetLabel();

    // use smearing matrix to smear the original parameters
    SmearTrack(xx,xxsm,CovMatSmear);

    AliTPCtrack* tpctrack = new AliTPCtrack(0,xxsm,cc,xref,alpha);
    tpctrack->SetLabel(lab);

    // fill the array
    newtarray.AddLast(tpctrack);
  
  }

  DBfile->Close();

  return; 
}
//-----------------------------------------------------------------
Commit	Line	Data
6eb67451	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	$Log$
	17	*/
	18	/*************************************************************************
	19	* *
	20	* This class builds AliTPCtrack objects from generated tracks to feed *
	21	* ITS tracking (V2). The AliTPCtrack is built from its first hit in *
	22	* the TPC. The track is assigned a Kalman-like covariance matrix *
	23	* depending on its pT and pseudorapidity and track parameters are *
	24	* smeared according this covariance matrix. *
	25	* Output file contains sorted tracks, ready for matching with ITS *
	26	* *
	27	* Test macro is: AliBarrelRec_TPCparam.C *
	28	* *
	29	* Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it *
	30	* *
	31	**************************************************************************/
	32	#include "AliTPCtrackerParam.h"
	33	#include "alles.h"
	34	#include "AliMagF.h"
	35	#include "AliTPCtrack.h"
	36	#include "TMatrixD.h"
	37	#include "AliKalmanTrack.h"
	38
	39	ClassImp(AliTPCtrackerParam)
	40
	41	//-----------------------------------------------------------------
	42	AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz)
	43	{
	44	fColl = coll; // collision code (0: PbPb6000)
	45	fBz = Bz; // value of the z component of L3 field (Tesla)
	46
	47	}
	48	//-----------------------------------------------------------------
	49	AliTPCtrackerParam::~AliTPCtrackerParam()
	50	{}
	51	//-----------------------------------------------------------------
	52
	53	Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile inp, TFile out, Int_t n)
	54	{
	55
	56	if(fColl!=0) {
	57	cerr<<"AliTPCtrackerParam::BuildTPCtracks: Invalid collision!\n
	58	Available: 0 -> PbPb6000"<<endl; return 0;
	59	}
	60	if(fBz!=0.4) {
	61	cerr<<"AliTPCtrackerParam::BuildTPCtracks: Invalid field!\n
	62	Available: 0.4"<<endl; return 0;
	63	}
	64
65	TFile infile=(TFile)inp;
66
67	AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
68
69	// Get gAlice object from file
70	if(!(gAlice=(AliRun*)infile->Get("gAlice"))) {
71	cerr<<"gAlice has not been found on galice.root !\n";
72	return 1;
73	}
74
75	// Set random number generator seed
76	TDatime t;
77	gRandom->SetSeed(t.Get());
78
79
80	// loop over first n events in file
81	for(Int_t evt=0; evt<n; evt++){
82	cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
83
84	AliTPCtrack *tpctrack=0;
85
86	// tree for TPC tracks
87	Char_t tname[100];
88	sprintf(tname,"TreeT_TPC_%d",evt);
89	TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
90	tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);
91
92	// array for TPC tracks
93	TObjArray tarray(20000);
94
95
96	Int_t nparticles=gAlice->GetEvent(evt);
97
98	Bool_t done[500000];
99	for(Int_t l=0; l<500000; l++) { done[l]=kFALSE; }
100
101
102	// Get TPC detector
103	AliTPC TPC=(AliTPC)gAlice->GetDetector("TPC");
104	Int_t ver = TPC->IsVersion();
105	cerr<<"+++ TPC version "<<ver<<" has been found !\n";
106	AliTPCParam digp=(AliTPCParam)infile->Get("75x40_100x60");
107	if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
108	TPC->SetParam(digp);
109
110	// Get TreeH with hits
111	TTree *TH=gAlice->TreeH();
112	Int_t ntracks=(Int_t)TH->GetEntries();
113	cerr<<"+++\n+++ Number of particles in event "<<evt<<": "<<nparticles<<"\n+++\n+++ Number of \"primary tracks\" (entries in TreeH): "<<ntracks<<"\n+++\n\n";
114
115	TParticle* Part;
116	AliTPChit* tpcHit;
117	Double_t hPx,hPy,hPz,hPt,xg,yg,zg,xl,yl,zl;
118	Double_t Alpha;
119	Float_t cosAlpha,sinAlpha;
120	Int_t label,Pdg,charge,bin;
121	Int_t tracks=0;
122	//Int_t N1=0,N2=0;
123
124	// loop over entries in TreeH
125	for(Int_t i=0; i<ntracks; i++) {
126	if(i%1000==0) cerr<<" --- Processing primary track "<<i<<" of "<<ntracks<<" ---\r";
127	TPC->ResetHits();
128	TH->GetEvent(i);
129	// Get FirstHit
130	tpcHit=(AliTPChit*)TPC->FirstHit(-1);
131	for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) {
132	if(tpcHit->fQ !=0.) continue;
133	// Get particle momentum at hit
134	hPx=tpcHit->X(); hPy=tpcHit->Y(); hPz=tpcHit->Z();
135	hPt=TMath::Sqrt(hPxhPx+hPyhPy);
136	// reject hits with Pt<mag*0.45 GeV/c
137	if(hPt<(fBz*0.45)) continue;
138
139	// Get track label
140	label=tpcHit->Track();
141	// check if this track has already been processed
142	if(done[label]) continue;
143	// electric charge
144	Part = gAlice->Particle(label);
145	Pdg = Part->GetPdgCode();
146	if(Pdg>200 \|\| Pdg==-11 \|\| Pdg==-13) { charge=1; }
147	else if(Pdg<-200 \|\| Pdg==11 \|\| Pdg==13) { charge=-1; }
148	else continue;
149
150
151	if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break;
152	if(tpcHit->fQ != 0.) continue;
153	// Get global coordinates of hit
154	xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
155	if(TMath::Sqrt(xgxg+ygyg)>90.) continue;
156
157	// Get TPC sector, Alpha angle and local coordinates
158	// printf("Sector %d\n",tpcHit->fSector);
159	digp->AdjustCosSin(tpcHit->fSector,cosAlpha,sinAlpha);
160	Alpha = TMath::ATan2(sinAlpha,cosAlpha);
161	xl = xgcosAlpha + ygsinAlpha;
162	yl =-xgsinAlpha + ygcosAlpha;
163	zl = zg;
164	//printf("Alpha %f xl %f yl %f zl %f\n",Alpha,xl,yl,zl);
165
166	// reject tracks which are not in the TPC acceptance
167	if(TMath::Abs(zl+(244.-xl)*hPz/hPt)>252.) continue;
168
169	// Get bin in pT,eta
170	bin = GetBin(hPt,Part->Eta());
171
172	// Apply selection according to TPC efficiency
173	//if(abs(Pdg)==211) N1++;
174	if(!SelectedTrack(Pdg,hPt,Part->Eta())) continue;
175	//if(abs(Pdg)==211) N2++;
176
177	// Mark track as "done"
178	done[label]=kTRUE; tracks++;
179
180	// create AliTPCtrack object
181	tpctrack = BuildTrack(Alpha,xl,yl,zl,hPx,hPy,hPz,hPt,charge,label);
182
183	// put track in array
184	tarray.AddLast(tpctrack);
185
186	}
187
188	} // loop over entries in TreeH
189
190	TObjArray newtarray(20000);
191
192	// assing covariance matrixes and smear track parameters
193	CookTracks(tarray,newtarray);
194
195	// sort array with TPC tracks (decreasing pT)
196	newtarray.Sort();
197
198
199	Int_t arrentr = newtarray.GetEntriesFast();
200	//printf("\n %d \n\n",arrentr);
201	for(Int_t l=0; l<arrentr; l++) {
202	tpctrack=(AliTPCtrack*)newtarray.UncheckedAt(l);
203	tracktree->Fill();
204	}
205
206	// write the tree with tracks in the output file
207	out->cd();
208	tracktree->Write();
209
210	delete tracktree;
211
212	printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
213	//printf("Average Eff: %f\n",(Float_t)N2/N1);
214
215	} // loop on events
216
217	return 0;
218	}
219
220	//-----------------------------------------------------------------
221	AliTPCtrack* AliTPCtrackerParam::BuildTrack(Double_t alpha,Double_t x,
222	Double_t y,Double_t z,Double_t Px,
223	Double_t Py,Double_t Pz,Double_t Pt,
224	Int_t ch,Int_t lab)
225	{
226	Double_t xref = x;
227	Double_t xx[5],cc[15];
228	cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.;
229	cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
230
231	// Magnetic field
232	TVector3 Bfield(0.,0.,fBz);
233
234
235	// radius [cm] of track projection in (x,y)
236	Double_t rho = Pt*100/0.299792458/Bfield.Z();
237	// center of track projection in local reference frame
238	TVector3 MomH,PosH;
239
240
241	// position (local) and momentum (local) at the hit
242	// in the bending plane (z=0)
243	PosH.SetXYZ(x,y,0.);
244	MomH.SetXYZ(PxTMath::Cos(alpha)+PyTMath::Sin(alpha),-PxTMath::Sin(alpha)+PyTMath::Cos(alpha),0.);
245	TVector3 Vrho = MomH.Cross(Bfield);
246	Vrho *= ch;
247	Vrho.SetMag(rho);
248
249	TVector3 Vcenter = PosH+Vrho;
250
251	Double_t x0 = Vcenter.X();
252
253	// fX = xref X-coordinate of this track (reference plane)
254	// fAlpha = Alpha Rotation angle the local (TPC sector)
255	// fP0 = YL Y-coordinate of a track
256	// fP1 = ZG Z-coordinate of a track
257	// fP2 = C*x0 x0 is center x in rotated frame
258	// fP3 = Tgl tangent of the track momentum dip angle
259	// fP4 = C track curvature
260	xx[0] = y;
261	xx[1] = z;
262	xx[3] = Pz/Pt;
263	xx[4] = -ch/rho;
264	xx[2] = xx[4]*x0;
265
266	// create the object AliTPCtrack
267	AliTPCtrack* track = new AliTPCtrack(0,xx,cc,xref,alpha);
268	// set the label
269	track->SetLabel(lab);
270
271	return track;
272	}
273
274	//-----------------------------------------------------------------
275	Bool_t AliTPCtrackerParam::SelectedTrack(Int_t pdg,Double_t pt,Double_t eta)
276	{
277	Double_t eff=0.;
278
279	//eff computed with \| zl+(244-xl)*pz/pt \| < 252
280	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};
281	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};
282	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};
283	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};
284	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};
285
286
287	if(TMath::Abs(pdg)==211) eff = LinearInterpolation(9,effPi,pt,eta);
288	if(TMath::Abs(pdg)==321) eff = LinearInterpolation(6,effK,pt,eta);
289	if(TMath::Abs(pdg)==2212) eff = LinearInterpolation(5,effP,pt,eta);
290	if(TMath::Abs(pdg)==11) eff = LinearInterpolation(5,effEl,pt,eta);
291	if(TMath::Abs(pdg)==13) eff = LinearInterpolation(5,effMu,pt,eta);
292
293	if(gRandom->Rndm() < eff) return kTRUE;
294
295	return kFALSE;
296	}
297
298	//-----------------------------------------------------------------
299	Double_t AliTPCtrackerParam::LinearInterpolation(Int_t ptBins,Double_t *value,
300	Double_t pt,Double_t eta)
301	{
302	Double_t kValue=0,kValue1=0,kValue2=0;
303	Int_t kEtaSide = (TMath::Abs(eta)<.45 ? 0 : 1);
304	Double_t kEta[3]={0.15,0.45,0.75};
305	Double_t kPt[9]={0.244,0.390,0.676,1.190,2.36,4.,6.,10.,20.};
306	if(ptBins==6) kPt[5]=10.;
307
308	for(Int_t i=0; i<ptBins; i++) {
309	if(pt<kPt[i]) {
310	if(i==0) i=1;
311	kValue1 = value[3i-3+kEtaSide]+(value[3i-3+kEtaSide]-value[3i+kEtaSide])/(kPt[i-1]-kPt[i])(pt-kPt[i-1]);
312	kValue2 = value[3i-3+kEtaSide+1]+(value[3i-3+kEtaSide+1]-value[3i+kEtaSide+1])/(kPt[i-1]-kPt[i])(pt-kPt[i-1]);
313
314	kValue = kValue1+(kValue1-kValue2)/(kEta[kEtaSide]-kEta[kEtaSide+1])*(eta-kEta[kEtaSide]);
315	break;
316	}
317	if(i==ptBins-1) {
318	i=ptBins-2;
319	kValue1 = value[3i-3+kEtaSide]+(value[3i-3+kEtaSide]-value[3i+kEtaSide])/(kPt[i-1]-kPt[i])(pt-kPt[i-1]);
320	kValue2 = value[3i-3+kEtaSide+1]+(value[3i-3+kEtaSide+1]-value[3i+kEtaSide+1])/(kPt[i-1]-kPt[i])(pt-kPt[i-1]);
321
322	kValue = kValue1+(kValue1-kValue2)/(kEta[kEtaSide]-kEta[kEtaSide+1])*(eta-kEta[kEtaSide]);
323	break;
324	}
325	}
326
327	return kValue;
328	}
329
330	//-----------------------------------------------------------------
331	Int_t AliTPCtrackerParam::GetBin(Double_t pt,Double_t eta) {
332
333	if(TMath::Abs(eta)<0.3) {
334	if(pt<0.3) return 0;
335	if(pt>=0.3 && pt<0.5) return 3;
336	if(pt>=0.5 && pt<1.) return 6;
337	if(pt>=1. && pt<1.5) return 9;
338	if(pt>=1.5 && pt<3.) return 12;
339	if(pt>=3. && pt<5.) return 15;
340	if(pt>=5. && pt<7.) return 18;
341	if(pt>=7. && pt<15.) return 21;
342	if(pt>=15.) return 24;
343	}
344	if(TMath::Abs(eta)>=0.3 && TMath::Abs(eta)<0.6) {
345	if(pt<0.3) return 1;
346	if(pt>=0.3 && pt<0.5) return 4;
347	if(pt>=0.5 && pt<1.) return 7;
348	if(pt>=1. && pt<1.5) return 10;
349	if(pt>=1.5 && pt<3.) return 13;
350	if(pt>=3. && pt<5.) return 16;
351	if(pt>=5. && pt<7.) return 19;
352	if(pt>=7. && pt<15.) return 22;
353	if(pt>=15.) return 25;
354	}
355	if(TMath::Abs(eta)>=0.6) {
356	if(pt<0.3) return 2;
357	if(pt>=0.3 && pt<0.5) return 5;
358	if(pt>=0.5 && pt<1.) return 8;
359	if(pt>=1. && pt<1.5) return 11;
360	if(pt>=1.5 && pt<3.) return 14;
361	if(pt>=3. && pt<5.) return 17;
362	if(pt>=5. && pt<7.) return 20;
363	if(pt>=7. && pt<15.) return 23;
364	if(pt>=15.) return 26;
365	}
366
367	return -1;
368
369	}
370
371	//-----------------------------------------------------------------
372
373	TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,Double_t eta)
374	{
375	TMatrixD CovMat(5,5);
376
377	CovMat(0,0)=cc[0];
378	CovMat(1,0)=cc[1]; CovMat(0,1)=CovMat(1,0);
379	CovMat(1,1)=cc[2];
380	CovMat(2,0)=cc[3]; CovMat(0,2)=CovMat(2,0);
381	CovMat(2,1)=cc[4]; CovMat(1,2)=CovMat(2,1);
382	CovMat(2,2)=cc[5];
383	CovMat(3,0)=cc[6]; CovMat(0,3)=CovMat(3,0);
384	CovMat(3,1)=cc[7]; CovMat(1,3)=CovMat(3,1);
385	CovMat(3,2)=cc[8]; CovMat(2,3)=CovMat(3,2);
386	CovMat(3,3)=cc[9];
387	CovMat(4,0)=cc[10]; CovMat(0,4)=CovMat(4,0);
388	CovMat(4,1)=cc[11]; CovMat(1,4)=CovMat(4,1);
389	CovMat(4,2)=cc[12]; CovMat(2,4)=CovMat(4,2);
390	CovMat(4,3)=cc[13]; CovMat(3,4)=CovMat(4,3);
391	CovMat(4,4)=cc[14];
392
393	TMatrixD ScaleMat(5,5);
394	for(Int_t k=0;k<5;k++) {
395	for(Int_t l=0;l<5;l++) {
396	ScaleMat(k,l)=0.;
397	}
398	}
399
400
401	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};
402	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};
403	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};
404	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};
405	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};
406
407
408	ScaleMat(0,0) = LinearInterpolation(9,s00,pt,eta);
409	ScaleMat(1,1) = LinearInterpolation(9,s11,pt,eta);
410	ScaleMat(2,2) = LinearInterpolation(9,s22,pt,eta);
411	ScaleMat(3,3) = LinearInterpolation(9,s33,pt,eta);
412	ScaleMat(4,4) = LinearInterpolation(9,s44,pt,eta);
413
414	TMatrixD Mat(ScaleMat,TMatrixD::kMult,CovMat);
415	TMatrixD CovMatSmear(Mat,TMatrixD::kMult,ScaleMat);
416
417
418	return CovMatSmear;
419	}
420
421	//-----------------------------------------------------------------
422	void AliTPCtrackerParam::SmearTrack(Double_t* xx,Double_t* xxsm,TMatrixD cov) {
423
424	// build triangular matrix Bmat
425	TMatrixD* Bmat = new TMatrixD(5,5);
426
427	for(Int_t j=0;j<5;j++){
428	Double_t accum = 0;
429	for(Int_t k=0;k<j;k++){
430	accum += (Bmat)(j,k) (*Bmat)(j,k);
431	}
432	(*Bmat)(j,j)=TMath::Sqrt(TMath::Abs(cov(j,j)-accum));
433	for(Int_t i=j+1;i<5;i++){
434	accum = 0;
435	for(Int_t k=0;k<j;k++){
436	accum+=(Bmat)(i,k) (*Bmat)(j,k);
437	}
438	(Bmat)(i,j) = (cov(i,j)-accum) / (Bmat)(j,j);
439	}
440	}
441
442	// get array of numbers with normal distribution
443	TArrayD norm;
444	norm.Set(5);
445	for(Int_t l=0;l<5;l++) {
446	norm[l]=gRandom->Gaus();
447	}
448
449	// use Bmat matrix to generate correlated numbers
450	TArrayD corr;
451	corr.Set(5);
452	for(Int_t i=0;i<5;i++){
453	corr[i]=0;
454	for(Int_t j=0;j<=i;j++) corr[i] += (Bmat)(i,j)norm[j];
455	}
456
457	delete Bmat;
458
459	for(Int_t l=0;l<5;l++) {
460	xxsm[l] = xx[l]+corr[l];
461	}
462
463	return;
464	}
465
466	//-----------------------------------------------------------------
467	void AliTPCtrackerParam::CookTracks(TObjArray& tarray,TObjArray& newtarray)
468	{
469	TString* s = new TString("$ALICE_ROOT/TPC/CovMatrixDB_");
470	if(fColl==0) s->Append("PbPb6000");
471	if(fBz==0.4) s->Append("_B0.4T.root");
472	// open file with matrixes DB
473	TFile* DBfile = new TFile(s->Data());
474
475	AliTPCtrack* track=new AliTPCtrack;
476	Int_t entr = (Int_t)tarray.GetEntriesFast();
477
478	for(Int_t k=0; k<entr; k++) {
479	track=(AliTPCtrack*)tarray.UncheckedAt(k);
480
481	Double_t pt = 1/TMath::Abs(track->Get1Pt());
482	Double_t eta = -TMath::Log(TMath::Tan(0.25TMath::Pi()-0.5TMath::ATan(track->GetTgl())));
483
484	Int_t bin = GetBin(pt,eta);
485
486	// string with tree name
487	TString str("Tree_bin");
488	str+=bin;
489
490	// get the right tree
491	TTree* CovTree = (TTree*)DBfile->Get(str.Data());
492	TArrayF* matrix = new TArrayF;
493	CovTree->SetBranchAddress("matrixes",&matrix);
494
495	// get random entry from the tree
496	Int_t Entries = (Int_t)CovTree->GetEntries();
497	CovTree->GetEvent(gRandom->Integer(Entries));
498
499	Double_t xref,alpha,xx[5],xxsm[5],cc[15];
500	Int_t lab;
501	// get P and Cosl from track
502	Double_t Cosl = TMath::Cos(TMath::ATan(track->GetTgl()));
503	Double_t P = 1/TMath::Abs(track->Get1Pt())/Cosl;
504
505	// get covariance matrix from regularized matrix
506	cc[0]=matrix->At(0)*(1.588e-3+1.911e-4/TMath::Power(P,1.5));
507	cc[1]=matrix->At(1);
508	cc[2]=matrix->At(2)*(1.195e-3+0.8102e-3/P);
509	cc[3]=matrix->At(3)*(7.275e-5+1.181e-5/TMath::Power(P,1.5));
510	cc[4]=matrix->At(4);
511	cc[5]=matrix->At(5)*(5.163e-6+1.138e-6/TMath::Power(P,2.)/Cosl);
512	cc[6]=matrix->At(6);
513	cc[7]=matrix->At(7)*(1.176e-5+1.175e-5/TMath::Power(P,1.5)/Cosl/Cosl/Cosl);
514	cc[8]=matrix->At(8);
515	cc[9]=matrix->At(9)*(1.289e-7+4.636e-7/TMath::Power(P,1.7)/Cosl/Cosl/Cosl/Cosl);
516	cc[10]=matrix->At(10)*(4.056e-7+4.379e-8/TMath::Power(P,1.5));
517	cc[11]=matrix->At(11);
518	cc[12]=matrix->At(12)*(3.049e-8+8.244e-9/TMath::Power(P,2.)/Cosl);
519	cc[13]=matrix->At(13);
520	cc[14]=matrix->At(14)*(1.847e-10+5.822e-11/TMath::Power(P,2.)/Cosl/Cosl/Cosl);
521
522	TMatrixD CovMatSmear(5,5);
523
524	CovMatSmear = GetSmearingMatrix(cc,pt,eta);
525
526	// get track original parameters
527	xref=track->GetX();
528	alpha=track->GetAlpha();
529	xx[0]=track->GetY();
530	xx[1]=track->GetZ();
531	xx[2]=track->GetX()*track->GetC()-track->GetSnp();
532	xx[3]=track->GetTgl();
533	xx[4]=track->GetC();
534	lab=track->GetLabel();
535
536	// use smearing matrix to smear the original parameters
537	SmearTrack(xx,xxsm,CovMatSmear);
538
539	AliTPCtrack* tpctrack = new AliTPCtrack(0,xxsm,cc,xref,alpha);
540	tpctrack->SetLabel(lab);
541
542	// fill the array
543	newtarray.AddLast(tpctrack);
544
545	}
546
547	DBfile->Close();
548
549	return;
550	}
551	//-----------------------------------------------------------------
552
553
554
555
556
557
558
559
560
561
562
563