New class by Andrea Dainese. It deals with the parameters used by
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerParam.cxx
CommitLineData
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. *
c68184b6 14 **************************************************************************/
15
16/*
17$Log$
18*/
19
20
21 /**************************************************************************
6eb67451 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 *
70521312 27 * smeared according to this covariance matrix. *
6eb67451 28 * Output file contains sorted tracks, ready for matching with ITS *
29 * *
e130146c 30 * For details: *
31 * http://www.pd.infn.it/alipd/talks/soft/adIII02/TPCtrackingParam.htm *
32 * *
33 * Test macro is: AliBarrelRec_TPCparam.C *
6eb67451 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"
70521312 44#include "AliMagFCM.h"
45#include "AliGausCorr.h"
46
47
6eb67451 48ClassImp(AliTPCtrackerParam)
49
50//-----------------------------------------------------------------
51AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz)
52{
e130146c 53//-----------------------------------------------------------------
54// This is the class conctructor
55//-----------------------------------------------------------------
56
6eb67451 57 fColl = coll; // collision code (0: PbPb6000)
58 fBz = Bz; // value of the z component of L3 field (Tesla)
59
60}
61//-----------------------------------------------------------------
62AliTPCtrackerParam::~AliTPCtrackerParam()
63{}
64//-----------------------------------------------------------------
65
66Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out, Int_t n)
67{
e130146c 68//-----------------------------------------------------------------
69// This function creates the TPC parameterized tracks
70//-----------------------------------------------------------------
71
6eb67451 72 if(fColl!=0) {
f3d60f45 73 cerr<<"AliTPCtrackerParam::BuildTPCtracks: Invalid collision!\n";
74 cerr<<" Available: 0 -> PbPb6000"<<endl; return 0;
6eb67451 75 }
76 if(fBz!=0.4) {
f3d60f45 77 cerr<<"AliTPCtrackerParam::BuildTPCtracks: Invalid field!\n";
78 cerr<<" Available: 0.4"<<endl; return 0;
6eb67451 79 }
80
81 TFile *infile=(TFile*)inp;
82
6eb67451 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
70521312 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);
6eb67451 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");
7a09f434 129 if(digp){
130 delete digp;
131 digp = new AliTPCParamSR();
132 }
133 else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60");
134
6eb67451 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;
e130146c 146 Double_t alpha;
6eb67451 147 Float_t cosAlpha,sinAlpha;
e130146c 148 Int_t label,pdg,charge,bin;
6eb67451 149 Int_t tracks=0;
e130146c 150 //Int_t nSel=0,nAcc=0;
6eb67451 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);
e130146c 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; }
6eb67451 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);
e130146c 188 alpha = TMath::ATan2(sinAlpha,cosAlpha);
6eb67451 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
e130146c 201 //if(TMath::Abs(pdg)==211) nAcc++;
202 if(!SelectedTrack(pdg,hPt,Part->Eta())) continue;
203 //if(TMath::Abs(pdg)==211) nSel++;
6eb67451 204
205 // Mark track as "done"
206 done[label]=kTRUE; tracks++;
207
208 // create AliTPCtrack object
e130146c 209 tpctrack = BuildTrack(alpha,xl,yl,zl,hPx,hPy,hPz,hPt,charge,label);
6eb67451 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);
e130146c 241 //printf("Average Eff: %f\n",(Float_t)nSel/nAcc);
6eb67451 242
243 } // loop on events
244
245 return 0;
246}
247
248//-----------------------------------------------------------------
249AliTPCtrack* AliTPCtrackerParam::BuildTrack(Double_t alpha,Double_t x,
e130146c 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
6eb67451 253{
e130146c 254//-----------------------------------------------------------------
255// This function uses GEANT info to set true track parameters
256//-----------------------------------------------------------------
6eb67451 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
e130146c 263 TVector3 bfield(0.,0.,fBz);
6eb67451 264
265
266 // radius [cm] of track projection in (x,y)
e130146c 267 Double_t rho = pt*100./0.299792458/bfield.Z();
6eb67451 268 // center of track projection in local reference frame
e130146c 269 TVector3 hmom,hpos;
6eb67451 270
271
272 // position (local) and momentum (local) at the hit
273 // in the bending plane (z=0)
e130146c 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);
6eb67451 279
e130146c 280 TVector3 vcenter = hpos+vrho;
6eb67451 281
e130146c 282 Double_t x0 = vcenter.X();
6eb67451 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;
e130146c 293 xx[3] = pz/pt;
6eb67451 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//-----------------------------------------------------------------
306Bool_t AliTPCtrackerParam::SelectedTrack(Int_t pdg,Double_t pt,Double_t eta)
e130146c 307 const {
308//-----------------------------------------------------------------
309// This function makes a selection according to TPC tracking efficiency
310//-----------------------------------------------------------------
311
6eb67451 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//-----------------------------------------------------------------
334Double_t AliTPCtrackerParam::LinearInterpolation(Int_t ptBins,Double_t *value,
e130146c 335 Double_t trkPt,Double_t trkEta) const
6eb67451 336{
e130146c 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.;
6eb67451 345
346 for(Int_t i=0; i<ptBins; i++) {
e130146c 347 if(trkPt<pt[i]) {
6eb67451 348 if(i==0) i=1;
e130146c 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]);
6eb67451 351
e130146c 352 intValue = intValue1+(intValue1-intValue2)/(eta[etaSide]-eta[etaSide+1])*(trkEta-eta[etaSide]);
6eb67451 353 break;
354 }
355 if(i==ptBins-1) {
356 i=ptBins-2;
e130146c 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]);
6eb67451 359
e130146c 360 intValue = intValue1+(intValue1-intValue2)/(eta[etaSide]-eta[etaSide+1])*(trkEta-eta[etaSide]);
6eb67451 361 break;
362 }
363 }
364
e130146c 365 return intValue;
6eb67451 366}
367
368//-----------------------------------------------------------------
e130146c 369Int_t AliTPCtrackerParam::GetBin(Double_t pt,Double_t eta) const {
370//-----------------------------------------------------------------
371// This function tells bin number in a grid (pT,eta)
372//-----------------------------------------------------------------
6eb67451 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//-----------------------------------------------------------------
e130146c 412TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,Double_t eta) const
6eb67451 413{
e130146c 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);
6eb67451 436 for(Int_t k=0;k<5;k++) {
437 for(Int_t l=0;l<5;l++) {
e130146c 438 stretchMat(k,l)=0.;
6eb67451 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
e130146c 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);
6eb67451 455
e130146c 456 TMatrixD mat(stretchMat,TMatrixD::kMult,covMat);
457 TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat);
6eb67451 458
459
e130146c 460 return covMatSmear;
6eb67451 461}
462
463//-----------------------------------------------------------------
e130146c 464void AliTPCtrackerParam::SmearTrack(Double_t* xx,Double_t* xxsm,TMatrixD cov) const {
465//-----------------------------------------------------------------
466// This function smears track parameters according to streched cov. matrix
467//-----------------------------------------------------------------
70521312 468 AliGausCorr *corgen = new AliGausCorr(cov,5);
469 TArrayD corr(5);
470 corgen->GetGaussN(corr);
471 delete corgen;
472 corgen = 0;
6eb67451 473
474 for(Int_t l=0;l<5;l++) {
475 xxsm[l] = xx[l]+corr[l];
476 }
477
478 return;
479}
480
481//-----------------------------------------------------------------
482void AliTPCtrackerParam::CookTracks(TObjArray& tarray,TObjArray& newtarray)
e130146c 483const {
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");
6eb67451 491 // open file with matrixes DB
e130146c 492 TFile* fileDB = TFile::Open(s.Data());
6eb67451 493
70521312 494 AliTPCtrack* track = 0;
e130146c 495 Int_t arrayEntries = (Int_t)tarray.GetEntriesFast();
6eb67451 496
e130146c 497 for(Int_t k=0; k<arrayEntries; k++) {
6eb67451 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
e130146c 510 TTree* covTree = (TTree*)fileDB->Get(str.Data());
6eb67451 511 TArrayF* matrix = new TArrayF;
e130146c 512 covTree->SetBranchAddress("matrixes",&matrix);
6eb67451 513
514 // get random entry from the tree
e130146c 515 Int_t treeEntries = (Int_t)covTree->GetEntries();
516 covTree->GetEvent(gRandom->Integer(treeEntries));
6eb67451 517
518 Double_t xref,alpha,xx[5],xxsm[5],cc[15];
519 Int_t lab;
520 // get P and Cosl from track
e130146c 521 Double_t cosl = TMath::Cos(TMath::ATan(track->GetTgl()));
522 Double_t p = 1/TMath::Abs(track->Get1Pt())/cosl;
6eb67451 523
524 // get covariance matrix from regularized matrix
e130146c 525 cc[0]=matrix->At(0)*(1.588e-3+1.911e-4/TMath::Power(p,1.5));
6eb67451 526 cc[1]=matrix->At(1);
e130146c 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));
6eb67451 529 cc[4]=matrix->At(4);
e130146c 530 cc[5]=matrix->At(5)*(5.163e-6+1.138e-6/TMath::Power(p,2.)/cosl);
6eb67451 531 cc[6]=matrix->At(6);
e130146c 532 cc[7]=matrix->At(7)*(1.176e-5+1.175e-5/TMath::Power(p,1.5)/cosl/cosl/cosl);
6eb67451 533 cc[8]=matrix->At(8);
e130146c 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));
6eb67451 536 cc[11]=matrix->At(11);
e130146c 537 cc[12]=matrix->At(12)*(3.049e-8+8.244e-9/TMath::Power(p,2.)/cosl);
6eb67451 538 cc[13]=matrix->At(13);
e130146c 539 cc[14]=matrix->At(14)*(1.847e-10+5.822e-11/TMath::Power(p,2.)/cosl/cosl/cosl);
6eb67451 540
e130146c 541 TMatrixD covMatSmear(5,5);
6eb67451 542
e130146c 543 covMatSmear = GetSmearingMatrix(cc,pt,eta);
6eb67451 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
e130146c 556 SmearTrack(xx,xxsm,covMatSmear);
6eb67451 557
558 AliTPCtrack* tpctrack = new AliTPCtrack(0,xxsm,cc,xref,alpha);
559 tpctrack->SetLabel(lab);
560
561 // fill the array
562 newtarray.AddLast(tpctrack);
70521312 563
e130146c 564 delete matrix;
70521312 565
6eb67451 566 }
567
e130146c 568 fileDB->Close();
569 delete fileDB;
70521312 570
571
6eb67451 572 return;
573}
574//-----------------------------------------------------------------
575
576
577
578
579
580
581
582
583
584
585
586