remove info message
[u/mrichter/AliRoot.git] / TRD / TRDbase / AliTRDtrackV1.cxx
CommitLineData
93fc2389 1
d9950a5a 2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
17/* $Id$ */
18
e2a1c98b 19#include "TVectorT.h"
203967fc 20#include "AliLog.h"
eb38ed55 21#include "AliESDtrack.h"
41880fac 22#include "AliTracker.h"
eb38ed55 23
d9950a5a 24#include "AliTRDtrackV1.h"
25#include "AliTRDcluster.h"
26#include "AliTRDcalibDB.h"
3afdab72 27#include "AliTRDReconstructor.h"
9dcc64cc 28#include "AliTRDPIDResponse.h"
d9950a5a 29#include "AliTRDrecoParam.h"
6951a056 30#include "AliTRDdEdxBaseUtils.h"
31#include "AliTRDdEdxCalibHistArray.h"
32#include "AliTRDdEdxCalibUtils.h"
33#include "AliTRDdEdxReconUtils.h"
d9950a5a 34
d9950a5a 35ClassImp(AliTRDtrackV1)
36
0906e73e 37///////////////////////////////////////////////////////////////////////////////
38// //
39// Represents a reconstructed TRD track //
40// Local TRD Kalman track //
41// //
42// Authors: //
43// Alex Bercuci <A.Bercuci@gsi.de> //
44// Markus Fasel <M.Fasel@gsi.de> //
45// //
46///////////////////////////////////////////////////////////////////////////////
d9950a5a 47
48//_______________________________________________________________
3b57a3f7 49AliTRDtrackV1::AliTRDtrackV1() : AliKalmanTrack()
eb2b4f91 50 ,fStatus(0)
352cef8f 51 ,fESDid(0)
3b57a3f7 52 ,fDE(0.)
e2a1c98b 53 ,fTruncatedMean(0)
d0037ea4 54 ,fNchamberdEdx(0)
55 ,fNclusterdEdx(0)
1f97f376 56 ,fNdEdxSlices(0)
4d6aee34 57 ,fkReconstructor(NULL)
58 ,fBackupTrack(NULL)
59 ,fTrackLow(NULL)
60 ,fTrackHigh(NULL)
0906e73e 61{
d9950a5a 62 //
63 // Default constructor
64 //
6e49cfdb 65 //printf("AliTRDtrackV1::AliTRDtrackV1()\n");
e44586fb 66
3b57a3f7 67 for(int i =0; i<3; i++) fBudget[i] = 0.;
68
69 Float_t pid = 1./AliPID::kSPECIES;
70 for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
71
72 for(int ip=0; ip<kNplane; ip++){
17896e82 73 fTrackletIndex[ip] = -1;
4d6aee34 74 fTracklet[ip] = NULL;
3b57a3f7 75 }
1f97f376 76 SetLabel(-123456789); // reset label
d9950a5a 77}
78
79//_______________________________________________________________
3b57a3f7 80AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) : AliKalmanTrack(ref)
eb2b4f91 81 ,fStatus(ref.fStatus)
352cef8f 82 ,fESDid(ref.fESDid)
3b57a3f7 83 ,fDE(ref.fDE)
e2a1c98b 84 ,fTruncatedMean(ref.fTruncatedMean)
d0037ea4 85 ,fNchamberdEdx(ref.fNchamberdEdx)
86 ,fNclusterdEdx(ref.fNclusterdEdx)
1f97f376 87 ,fNdEdxSlices(ref.fNdEdxSlices)
4d6aee34 88 ,fkReconstructor(ref.fkReconstructor)
89 ,fBackupTrack(NULL)
90 ,fTrackLow(NULL)
91 ,fTrackHigh(NULL)
d9950a5a 92{
93 //
3b57a3f7 94 // Copy constructor
d9950a5a 95 //
96
6e49cfdb 97 //printf("AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &)\n");
29b87567 98 SetBit(kOwner, kFALSE);
3b57a3f7 99 for(int ip=0; ip<kNplane; ip++){
100 fTrackletIndex[ip] = ref.fTrackletIndex[ip];
101 fTracklet[ip] = ref.fTracklet[ip];
102 }
e3cf3d02 103 if(ref.fTrackLow) fTrackLow = new AliExternalTrackParam(*ref.fTrackLow);
104 if(ref.fTrackHigh) fTrackHigh = new AliExternalTrackParam(*ref.fTrackHigh);
105
3b57a3f7 106 for (Int_t i = 0; i < 3;i++) fBudget[i] = ref.fBudget[i];
107
108 for(Int_t is = 0; is<AliPID::kSPECIES; is++) fPID[is] = ref.fPID[is];
109
110 AliKalmanTrack::SetNumberOfClusters(ref.GetNumberOfClusters());
d9950a5a 111}
112
113//_______________________________________________________________
3b57a3f7 114AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) : AliKalmanTrack()
eb2b4f91 115 ,fStatus(0)
352cef8f 116 ,fESDid(0)
3b57a3f7 117 ,fDE(0.)
e2a1c98b 118 ,fTruncatedMean(0)
d0037ea4 119 ,fNchamberdEdx(0)
120 ,fNclusterdEdx(0)
1f97f376 121 ,fNdEdxSlices(0)
4d6aee34 122 ,fkReconstructor(NULL)
123 ,fBackupTrack(NULL)
124 ,fTrackLow(NULL)
125 ,fTrackHigh(NULL)
d9950a5a 126{
127 //
3b57a3f7 128 // Constructor from AliESDtrack
d9950a5a 129 //
130
352cef8f 131 SetESDid(t.GetID());
3b57a3f7 132 SetLabel(t.GetLabel());
133 SetChi2(0.0);
0b433f72 134
1d26da6d 135 SetMass(t.GetMassForTracking());
3b57a3f7 136 AliKalmanTrack::SetNumberOfClusters(t.GetTRDncls());
17896e82 137 Int_t ti[]={-1, -1, -1, -1, -1, -1}; t.GetTRDtracklets(&ti[0]);
3b57a3f7 138 for(int ip=0; ip<kNplane; ip++){
17896e82 139 fTrackletIndex[ip] = ti[ip];
4d6aee34 140 fTracklet[ip] = NULL;
3b57a3f7 141 }
142 for(int i =0; i<3; i++) fBudget[i] = 0.;
143
144 Float_t pid = 1./AliPID::kSPECIES;
145 for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
146
147 const AliExternalTrackParam *par = &t;
4e6eb3f6 148 /* RS
88987d18 149 if (t.GetStatus() & AliESDtrack::kTRDbackup) {
3b57a3f7 150 par = t.GetOuterParam();
151 if (!par) {
152 AliError("No backup info!");
153 par = &t;
154 }
88987d18 155 }
4e6eb3f6 156 */
3b57a3f7 157 Set(par->GetX()
158 ,par->GetAlpha()
159 ,par->GetParameter()
160 ,par->GetCovariance());
161
162 if(t.GetStatus() & AliESDtrack::kTIME) {
163 StartTimeIntegral();
fc9b31a7 164 Double_t times[AliPID::kSPECIESC];
165 t.GetIntegratedTimes(times,AliPID::kSPECIESC);
3b57a3f7 166 SetIntegratedTimes(times);
167 SetIntegratedLength(t.GetIntegratedLength());
168 }
d9950a5a 169}
170
171//_______________________________________________________________
e17f4785 172AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 * const trklts, const Double_t p[5], const Double_t cov[15]
3b57a3f7 173 , Double_t x, Double_t alpha) : AliKalmanTrack()
eb2b4f91 174 ,fStatus(0)
352cef8f 175 ,fESDid(0)
3b57a3f7 176 ,fDE(0.)
e2a1c98b 177 ,fTruncatedMean(0)
d0037ea4 178 ,fNchamberdEdx(0)
179 ,fNclusterdEdx(0)
1f97f376 180 ,fNdEdxSlices(0)
4d6aee34 181 ,fkReconstructor(NULL)
182 ,fBackupTrack(NULL)
183 ,fTrackLow(NULL)
184 ,fTrackHigh(NULL)
d9950a5a 185{
0906e73e 186 //
187 // The stand alone tracking constructor
188 // TEMPORARY !!!!!!!!!!!
189 // to check :
190 // 1. covariance matrix
191 // 2. dQdl calculation
192 //
d9950a5a 193
68f9b6bd 194 Double_t b(GetBz());
195 Double_t cnv = (TMath::Abs(b) < 1.e-5) ? 1.e5 : 1./GetBz()/kB2C;
196
d9950a5a 197 Double_t pp[5] = { p[0]
68f9b6bd 198 , p[1]
199 , p[2]
200 , p[3]
201 , p[4]*cnv };
202
d9950a5a 203 Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5];
204 Double_t c32 = x*cov[13] - cov[ 8];
205 Double_t c20 = x*cov[10] - cov[ 3];
206 Double_t c21 = x*cov[11] - cov[ 4];
207 Double_t c42 = x*cov[14] - cov[12];
68f9b6bd 208
d9950a5a 209 Double_t cc[15] = { cov[ 0]
210 , cov[ 1], cov[ 2]
211 , c20, c21, c22
212 , cov[ 6], cov[ 7], c32, cov[ 9]
213 , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
68f9b6bd 214
ae3fbe1f 215 Double_t mostProbablePt=AliExternalTrackParam::GetMostProbablePt();
216 Double_t p0=TMath::Sign(1/mostProbablePt,pp[4]);
217 Double_t w0=cc[14]/(cc[14] + p0*p0), w1=p0*p0/(cc[14] + p0*p0);
68f9b6bd 218 AliDebug(3, Form("Pt mixing : w0[%4.2f] pt0[%5.3f] w1[%4.2f] pt[%5.3f]", w0, 1./p0, w1, 1./pp[4]));
219
ae3fbe1f 220 pp[4] = w0*p0 + w1*pp[4];
221 cc[10]*=w1; cc[11]*=w1; cc[12]*=w1; cc[13]*=w1; cc[14]*=w1;
68f9b6bd 222 Set(x,alpha,pp,cc);
0217fcd0 223 AliDebug(2, Form("Init @ x[%6.2f] pt[%5.3f]", x, 1./pp[4]));
41702fec 224 Int_t ncls = 0;
68f9b6bd 225 for(int iplane=0; iplane<kNplane; iplane++){
17896e82 226 fTrackletIndex[iplane] = -1;
68f9b6bd 227 if(!trklts[iplane].IsOK()) fTracklet[iplane] = NULL;
41702fec 228 else{
e17f4785 229 fTracklet[iplane] = &trklts[iplane];
41702fec 230 ncls += fTracklet[iplane]->GetN();
231 }
68f9b6bd 232 }
41702fec 233 AliKalmanTrack::SetNumberOfClusters(ncls);
234 for(int i =0; i<3; i++) fBudget[i] = 0.;
235
236 Float_t pid = 1./AliPID::kSPECIES;
237 for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
1f97f376 238 SetLabel(-123456789); // reset label
d9950a5a 239}
240
241//_______________________________________________________________
3b57a3f7 242AliTRDtrackV1::~AliTRDtrackV1()
bb56afff 243{
4d8f1bd9 244 // Clean up all objects allocated by the track during its lifetime.
245 AliDebug(2, Form("Deleting track[%d]\n fBackupTrack[%p] fTrackLow[%p] fTrackHigh[%p] Owner[%c].", fESDid, (void*)fBackupTrack, (void*)fTrackLow, (void*)fTrackHigh, TestBit(kOwner)?'y':'n'));
76b60503 246
4d6aee34 247 if(fBackupTrack) delete fBackupTrack; fBackupTrack = NULL;
e3cf3d02 248
4d6aee34 249 if(fTrackLow) delete fTrackLow; fTrackLow = NULL;
250 if(fTrackHigh) delete fTrackHigh; fTrackHigh = NULL;
bb56afff 251
76b60503 252 for(Int_t ip=0; ip<kNplane; ip++){
253 if(TestBit(kOwner) && fTracklet[ip]) delete fTracklet[ip];
4d6aee34 254 fTracklet[ip] = NULL;
17896e82 255 fTrackletIndex[ip] = -1;
3b57a3f7 256 }
257}
258
259//_______________________________________________________________
e08a5492 260AliTRDtrackV1 &AliTRDtrackV1::operator=(const AliTRDtrackV1 &t)
261{
262 //
263 // Assignment operator
264 //
265
266 if (this != &t) {
5c5d503a 267 AliKalmanTrack::operator=(t);
e08a5492 268 ((AliTRDtrackV1 &) t).Copy(*this);
269 }
270
271 return *this;
272
273}
274
275//_____________________________________________________________________________
276void AliTRDtrackV1::Copy(TObject &t) const
277{
278 //
279 // Copy function
280 //
281
282 ((AliTRDtrackV1 &) t).fStatus = fStatus;
283 ((AliTRDtrackV1 &) t).fESDid = fESDid;
284 ((AliTRDtrackV1 &) t).fDE = fDE;
285 ((AliTRDtrackV1 &) t).fkReconstructor = fkReconstructor;
286 ((AliTRDtrackV1 &) t).fBackupTrack = 0x0;
287 ((AliTRDtrackV1 &) t).fTrackLow = 0x0;
288 ((AliTRDtrackV1 &) t).fTrackHigh = 0x0;
289
290 for(Int_t ip = 0; ip < kNplane; ip++) {
291 ((AliTRDtrackV1 &) t).fTrackletIndex[ip] = fTrackletIndex[ip];
292 ((AliTRDtrackV1 &) t).fTracklet[ip] = fTracklet[ip];
293 }
294 if (fTrackLow) {
295 ((AliTRDtrackV1 &) t).fTrackLow = new AliExternalTrackParam(*fTrackLow);
296 }
297 if (fTrackHigh){
298 ((AliTRDtrackV1 &) t).fTrackHigh = new AliExternalTrackParam(*fTrackHigh);
299 }
300
301 for (Int_t i = 0; i < 3; i++) {
302 ((AliTRDtrackV1 &) t).fBudget[i] = fBudget[i];
303 }
304 for (Int_t is = 0; is < AliPID::kSPECIES; is++) {
305 ((AliTRDtrackV1 &) t).fPID[is] = fPID[is];
306 }
307
308}
309
310//_______________________________________________________________
e73baf29 311Int_t AliTRDtrackV1::CookLabel(Float_t wrong, Int_t *labs, Float_t *freq)
3b57a3f7 312{
e73baf29 313// Set MC label for this track
314// On demand i.e. if arrays "labs" and "freq" are allocated by user returns :
315// nlabs = the no of distinct labels
316// labs = array of distinct labels in decreasing order of frequency
317// freq = frequency of each label in decreasing order
bb2db46c 318
e73baf29 319 Int_t ncl(0);
320 if(!(ncl = GetNumberOfClusters())) return 0;
321
322 Int_t s[2][kMAXCLUSTERSPERTRACK];
bb56afff 323 for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
e73baf29 324 s[0][i] = -1;
325 s[1][i] = 0;
bb56afff 326 }
e73baf29 327
328 Int_t label(-123456789), nlabels(0);
329 AliTRDcluster *c(NULL);
330 for (Int_t ip(0); ip < AliTRDgeometry::kNlayer; ip++) {
68f9b6bd 331 if(fTrackletIndex[ip]<0 || !fTracklet[ip]) continue;
e73baf29 332 for (Int_t ic(0); ic < AliTRDseedV1::kNclusters; ic++) {
3b57a3f7 333 if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
e73baf29 334 for (Int_t k(0); k < 3; k++) {
335 if ((label = c->GetLabel(k)) < 0) continue;
336 Int_t j(0);
337 while(j < kMAXCLUSTERSPERTRACK){
338 if(s[0][j]!=label && s[1][j]!=0){j++; continue;}
339 if(!s[1][j]) nlabels++;
340 s[0][j] = label; s[1][j]++;
341 break;
3b57a3f7 342 }
343 }
344 }
345 }
e73baf29 346 //printf(" Found %4d labels\n", nlabels);
347 Float_t prob(1.);
348 if(!nlabels){
349 AliError(Form("No MC labels found for track %d.", fESDid));
350 return 0;
351 } else if(nlabels==1) {
352 label = s[0][0];
353 if(labs && freq){labs[0]=label; freq[0]=1.;}
354 } else {
355 Int_t idx[kMAXCLUSTERSPERTRACK];
356 TMath::Sort(nlabels, s[1], idx);
357 label = s[0][idx[0]]; prob = s[1][idx[0]]/Float_t(ncl);
358 if(labs && freq){
359 for (Int_t i(0); i<nlabels; i++){
360 labs[i] = s[0][idx[i]];
361 freq[i] = s[1][idx[i]]/Float_t(ncl);
362 }
363 }
bb56afff 364 }
e73baf29 365 SetLabel((1.-prob > wrong)?-label:label);
366 return nlabels;
bb56afff 367}
368
369//_______________________________________________________________
d9950a5a 370Bool_t AliTRDtrackV1::CookPID()
371{
2a3191bb 372//
373// Cook the PID information for the track by delegating the omonim function of the tracklets.
374// Computes the number of tracklets used. The tracklet information are considered independent.
375// For the moment no global track measurement of PID is performed as for example to estimate
376// bremsstrahlung probability based on global chi2 of the track.
377//
378// The status bit AliESDtrack::kTRDpid is set during the call of AliTRDtrackV1::UpdateESDtrack().The PID performance of the
379//TRD for tracks with 6 tacklets is displayed below.
380//Begin_Html
381//<img src="TRD/trackPID.gif">
382//End_Html
383//
9dcc64cc 384 const AliTRDPIDResponse *pidResponse = AliTRDcalibDB::Instance()->GetPIDResponse(fkReconstructor->GetRecoParam()->GetPIDmethod());
385 if(!pidResponse){
386 AliError("PID Response not available");
eb2b4f91 387 return kFALSE;
3b57a3f7 388 }
1f97f376 389 Double_t dEdx[kNplane * (Int_t)AliTRDseedV1::kNdEdxSlices] = {0.};
390 Float_t trackletP[kNplane] = {0.};
391
392 fNdEdxSlices = pidResponse->GetNumberOfSlices();
9dcc64cc 393 for(Int_t iseed = 0; iseed < kNplane; iseed++){
394 if(!fTracklet[iseed]) continue;
395 trackletP[iseed] = fTracklet[iseed]->GetMomentum();
1f97f376 396// if(pidResponse->GetPIDmethod() == AliTRDPIDResponse::kLQ1D){
9dcc64cc 397 dEdx[iseed] = fTracklet[iseed]->GetdQdl();
1f97f376 398/* } else {
399 fTracklet[iseed]->CookdEdx(fNdEdxSlices);
9dcc64cc 400 const Float_t *trackletdEdx = fTracklet[iseed]->GetdEdx();
1f97f376 401 for(Int_t islice = 0; islice < fNdEdxSlices; islice++){
402 dEdx[iseed*fNdEdxSlices + islice] = trackletdEdx[islice];
9dcc64cc 403 }
1f97f376 404 }*/
9dcc64cc 405 }
1f97f376 406 pidResponse->GetResponse(fNdEdxSlices, dEdx, trackletP, fPID);
e2a1c98b 407
9a8b0e85 408 static Int_t nprint = 0;
409 if(!nprint){
6951a056 410 AliTRDdEdxBaseUtils::PrintControl();
9a8b0e85 411 nprint++;
412 }
413
e2a1c98b 414 //do truncated mean
6951a056 415 AliTRDdEdxCalibUtils::SetObjArray(AliTRDcalibDB::Instance()->GetPHQ());
416 const Double_t mag = AliTRDdEdxBaseUtils::IsExBOn() ? GetBz() : -1;
417 const Double_t charge = AliTRDdEdxBaseUtils::IsExBOn() ? Charge() : -1;
d0037ea4 418 fTruncatedMean = CookTruncatedMean(0, mag, charge, kTRUE, fNchamberdEdx, fNclusterdEdx);
e2a1c98b 419
eb2b4f91 420 return kTRUE;
421}
422
423//___________________________________________________________
424UChar_t AliTRDtrackV1::GetNumberOfTrackletsPID() const
425{
426// Retrieve number of tracklets used for PID calculation.
427
e20bef2b 428 UChar_t nPID = 0;
3b57a3f7 429 for(int ip=0; ip<kNplane; ip++){
68f9b6bd 430 if(fTrackletIndex[ip]<0 || !fTracklet[ip]) continue;
3b57a3f7 431 if(!fTracklet[ip]->IsOK()) continue;
3b57a3f7 432
e20bef2b 433 nPID++;
0906e73e 434 }
e20bef2b 435 return nPID;
eb2b4f91 436}
437
3b57a3f7 438//_______________________________________________________________
0349cc67 439AliTRDcluster* AliTRDtrackV1::GetCluster(Int_t id)
440{
4d6aee34 441 // Get the cluster at a certain position in the track
0349cc67 442 Int_t n = 0;
443 for(Int_t ip=0; ip<kNplane; ip++){
444 if(!fTracklet[ip]) continue;
445 if(n+fTracklet[ip]->GetN() <= id){
446 n+=fTracklet[ip]->GetN();
447 continue;
448 }
4d6aee34 449 AliTRDcluster *c = NULL;
8d2bec9e 450 for(Int_t ic=AliTRDseedV1::kNclusters; ic--;){
0349cc67 451 if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
452
453 if(n<id){n++; continue;}
454 return c;
455 }
456 }
4d6aee34 457 return NULL;
0349cc67 458}
459
460//_______________________________________________________________
3b57a3f7 461Int_t AliTRDtrackV1::GetClusterIndex(Int_t id) const
462{
4d6aee34 463 // Get the cluster index at a certain position in the track
3b57a3f7 464 Int_t n = 0;
465 for(Int_t ip=0; ip<kNplane; ip++){
466 if(!fTracklet[ip]) continue;
76b60503 467 if(n+fTracklet[ip]->GetN() <= id){
3b57a3f7 468 n+=fTracklet[ip]->GetN();
469 continue;
470 }
8d2bec9e 471 for(Int_t ic=AliTRDseedV1::kNclusters; ic--;){
87f70a90 472 if(!(fTracklet[ip]->GetClusters(ic))) continue;
76b60503 473 if(n<id){n++; continue;}
3b57a3f7 474 return fTracklet[ip]->GetIndexes(ic);
475 }
476 }
477 return -1;
0906e73e 478}
d9950a5a 479
480//_______________________________________________________________
b72f4eaf 481Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt, Double_t *cov) const
d9950a5a 482{
ed15ef4f 483// Compute chi2 between tracklet and track. The value is calculated at the radial position of the track
484// equal to the reference radial position of the tracklet (see AliTRDseedV1)
485//
486// The chi2 estimator is computed according to the following formula
487// BEGIN_LATEX
488// #chi^{2}=(X_{trklt}-X_{track})(C_{trklt}+C_{track})^{-1}(X_{trklt}-X_{track})^{T}
489// END_LATEX
490// where X=(y z), the position of the track/tracklet in the yz plane
491//
492
b72f4eaf 493 Double_t p[2] = { trklt->GetY(), trklt->GetZ()};
494 trklt->GetCovAt(trklt->GetX(), cov);
ed15ef4f 495 return AliExternalTrackParam::GetPredictedChi2(p, cov);
3b57a3f7 496}
d9950a5a 497
22a4ab0c 498//_______________________________________________________________
eb2b4f91 499Int_t AliTRDtrackV1::GetSector() const
500{
501 return Int_t(GetAlpha()/AliTRDgeometry::GetAlpha() + (GetAlpha()>0. ? 0 : AliTRDgeometry::kNsector));
502}
503
504//_______________________________________________________________
203967fc 505Bool_t AliTRDtrackV1::IsEqual(const TObject *o) const
506{
4d6aee34 507 // Checks whether two tracks are equal
203967fc 508 if (!o) return kFALSE;
509 const AliTRDtrackV1 *inTrack = dynamic_cast<const AliTRDtrackV1*>(o);
510 if (!inTrack) return kFALSE;
511
eb2b4f91 512 //if ( fPIDquality != inTrack->GetPIDquality() ) return kFALSE;
203967fc 513
d8069611 514 if(memcmp(fPID, inTrack->fPID, AliPID::kSPECIES*sizeof(Double32_t))) return kFALSE;
515 if(memcmp(fBudget, inTrack->fBudget, 3*sizeof(Double32_t))) return kFALSE;
516 if(memcmp(&fDE, &inTrack->fDE, sizeof(Double32_t))) return kFALSE;
517 if(memcmp(&fFakeRatio, &inTrack->fFakeRatio, sizeof(Double32_t))) return kFALSE;
518 if(memcmp(&fChi2, &inTrack->fChi2, sizeof(Double32_t))) return kFALSE;
519 if(memcmp(&fMass, &inTrack->fMass, sizeof(Double32_t))) return kFALSE;
520 if( fLab != inTrack->fLab ) return kFALSE;
521 if( fN != inTrack->fN ) return kFALSE;
522 Double32_t l(0.), in(0.);
523 l = GetIntegratedLength(); in = inTrack->GetIntegratedLength();
524 if(memcmp(&l, &in, sizeof(Double32_t))) return kFALSE;
525 l=GetX(); in=inTrack->GetX();
526 if(memcmp(&l, &in, sizeof(Double32_t))) return kFALSE;
527 l = GetAlpha(); in = inTrack->GetAlpha();
528 if(memcmp(&l, &in, sizeof(Double32_t))) return kFALSE;
529 if(memcmp(GetParameter(), inTrack->GetParameter(), 5*sizeof(Double32_t))) return kFALSE;
530 if(memcmp(GetCovariance(), inTrack->GetCovariance(), 15*sizeof(Double32_t))) return kFALSE;
203967fc 531
532 for (Int_t iTracklet = 0; iTracklet < kNplane; iTracklet++){
533 AliTRDseedV1 *curTracklet = fTracklet[iTracklet];
534 AliTRDseedV1 *inTracklet = inTrack->GetTracklet(iTracklet);
535 if (curTracklet && inTracklet){
536 if (! curTracklet->IsEqual(inTracklet) ) {
537 curTracklet->Print();
538 inTracklet->Print();
539 return kFALSE;
540 }
541 } else {
542 // if one tracklet exists, and corresponding
543 // in other track doesn't - return kFALSE
544 if(inTracklet || curTracklet) return kFALSE;
545 }
546 }
547
548 return kTRUE;
549}
550
551//_______________________________________________________________
22a4ab0c 552Bool_t AliTRDtrackV1::IsElectron() const
553{
4d6aee34 554 if(GetPID(0) > fkReconstructor->GetRecoParam()->GetPIDThreshold(GetP())) return kTRUE;
22a4ab0c 555 return kFALSE;
556}
557
3b57a3f7 558
559//_____________________________________________________________________________
b453ef55 560Int_t AliTRDtrackV1::MakeBackupTrack()
3b57a3f7 561{
b453ef55 562//
563// Creates a backup track
564// TO DO update quality check of the track.
565//
566
567 Float_t occupancy(0.); Int_t n(0), ncls(0);
568 for(Int_t il(AliTRDgeometry::kNlayer); il--;){
569 if(!fTracklet[il]) continue;
570 n++;
5c5d503a 571 occupancy+=fTracklet[il]->GetTBoccupancy()/AliTRDseedV1::kNtb;
b453ef55 572 ncls += fTracklet[il]->GetN();
573 }
574 if(!n) return -1;
575 occupancy/=n;
576
577 //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);
578
579 Int_t failedCutId(0);
580 if(GetChi2()/n > 5.0) failedCutId=1;
581 if(occupancy < 0.7) failedCutId=2;
582 //if(ratio1 > 0.6) &&
583 //if(ratio0+ratio1 > 1.5) &&
584 if(GetNCross() != 0) failedCutId=3;
585 if(TMath::Abs(GetSnp()) > 0.85) failedCutId=4;
586 if(ncls < 20) failedCutId=5;
587
588 if(failedCutId){
589 AliDebug(2, Form("\n"
590 "chi2/tracklet < 5.0 [%c] %5.2f\n"
591 "occupancy > 0.7 [%c] %4.2f\n"
592 "NCross == 0 [%c] %d\n"
593 "Abs(snp) < 0.85 [%c] %4.2f\n"
594 "NClusters > 20 [%c] %d"
595 ,(GetChi2()/n<5.0)?'y':'n', GetChi2()/n
596 ,(occupancy>0.7)?'y':'n', occupancy
597 ,(GetNCross()==0)?'y':'n', GetNCross()
598 ,(TMath::Abs(GetSnp())<0.85)?'y':'n', TMath::Abs(GetSnp())
599 ,(ncls>20)?'y':'n', ncls
600 ));
601 return failedCutId;
602 }
3b57a3f7 603
604 if(fBackupTrack) {
605 fBackupTrack->~AliTRDtrackV1();
606 new(fBackupTrack) AliTRDtrackV1((AliTRDtrackV1&)(*this));
b453ef55 607 return 0;
3b57a3f7 608 }
609 fBackupTrack = new AliTRDtrackV1((AliTRDtrackV1&)(*this));
b453ef55 610 return 0;
d9950a5a 611}
612
3b57a3f7 613//_____________________________________________________________________________
f8a9723f 614Int_t AliTRDtrackV1::GetProlongation(Double_t xk, Double_t &y, Double_t &z) const
3b57a3f7 615{
616 //
617 // Find a prolongation at given x
530fb4e2 618 // Return -1 if it does not exist
3b57a3f7 619 //
620
621 Double_t bz = GetBz();
530fb4e2 622 if (!AliExternalTrackParam::GetYAt(xk,bz,y)) return -1;
623 if (!AliExternalTrackParam::GetZAt(xk,bz,z)) return -1;
3b57a3f7 624
625 return 1;
626
627}
628
629//_____________________________________________________________________________
15a2657d 630Bool_t AliTRDtrackV1::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
0906e73e 631{
632 //
3b57a3f7 633 // Propagates this track to a reference plane defined by "xk" [cm]
634 // correcting for the mean crossed material.
0906e73e 635 //
3b57a3f7 636 // "xx0" - thickness/rad.length [units of the radiation length]
637 // "xrho" - thickness*density [g/cm^2]
638 //
0906e73e 639
952051c5 640 if (TMath::Abs(xk - GetX())<AliTRDReconstructor::GetEpsilon()*0.1) return kTRUE; // 10% of the tracker precision
dc5e390a 641
307aac71 642 Double_t xyz0[3] = {GetX(), GetY(), GetZ()}, // track position BEFORE propagation
dc5e390a 643 b[3]; // magnetic field
307aac71 644 GetBxByBz(b);
dc5e390a 645 if(!AliExternalTrackParam::PropagateToBxByBz(xk,b)) return kFALSE;
646
307aac71 647 // local track position AFTER propagation
648 Double_t xyz1[3] = {GetX(), GetY(), GetZ()};
2f4384e6 649// printf("x0[%6.2f] -> x1[%6.2f] dx[%6.2f] rho[%f]\n", xyz0[0], xyz1[0], xyz0[0]-xk, xrho/TMath::Abs(xyz0[0]-xk));
dc5e390a 650 if(xyz0[0] < xk) {
5c5d503a 651 xrho = -xrho;
3b57a3f7 652 if (IsStartedTimeIntegral()) {
dc5e390a 653 Double_t l2 = TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])
654 + (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1])
655 + (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
6790f7d7 656 Double_t crv = AliExternalTrackParam::GetC(b[2]);
3b57a3f7 657 if (TMath::Abs(l2*crv) > 0.0001) {
658 // Make correction for curvature if neccesary
dc5e390a 659 l2 = 0.5 * TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])
660 + (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]));
3b57a3f7 661 l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
dc5e390a 662 l2 = TMath::Sqrt(l2*l2 + (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
3b57a3f7 663 }
664 AddTimeStep(l2);
665 }
666 }
35315297 667 if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0, xrho, fMass)) return kFALSE;
fd40f855 668
3b57a3f7 669
670 {
671
672 // Energy losses
673 Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
1d26da6d 674 if (fMass<0) p2 *= 4; // q=2
35315297 675 Double_t beta2 = p2 / (p2 + fMass*fMass);
3b57a3f7 676 if ((beta2 < 1.0e-10) ||
677 ((5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0)) {
678 return kFALSE;
679 }
680
681 Double_t dE = 0.153e-3 / beta2
682 * (TMath::Log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
683 * xrho;
1d26da6d 684 if (fMass<0) dE *= 4; // q=2
3b57a3f7 685 fBudget[0] += xrho;
686
687 /*
688 // Suspicious part - think about it ?
689 Double_t kinE = TMath::Sqrt(p2);
690 if (dE > 0.8*kinE) dE = 0.8 * kinE; //
691 if (dE < 0) dE = 0.0; // Not valid region for Bethe bloch
692 */
693
694 fDE += dE;
695
696 /*
697 // Suspicious ! I.B.
698 Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE)); // Energy loss fluctuation
35315297 699 Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+fMass*fMass)/(p2*p2);
3b57a3f7 700 fCcc += sigmac2;
701 fCee += fX*fX * sigmac2;
702 */
703
704 }
705
706 return kTRUE;
0906e73e 707}
3b57a3f7 708
87a7fa94 709//_____________________________________________________________________________
3b57a3f7 710Int_t AliTRDtrackV1::PropagateToR(Double_t r,Double_t step)
87a7fa94 711{
712 //
3b57a3f7 713 // Propagate track to the radial position
714 // Rotation always connected to the last track position
87a7fa94 715 //
716
3b57a3f7 717 Double_t xyz0[3];
718 Double_t xyz1[3];
719 Double_t y;
720 Double_t z;
721
722 Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
723 // Direction +-
724 Double_t dir = (radius > r) ? -1.0 : 1.0;
725
726 for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
727
728 GetXYZ(xyz0);
729 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
730 Rotate(alpha,kTRUE);
731 GetXYZ(xyz0);
e20bef2b 732 if(GetProlongation(x,y,z)<0) return -1;
3b57a3f7 733 xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
734 xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
735 xyz1[2] = z;
736 Double_t param[7];
83dea92e 737 if(AliTracker::MeanMaterialBudget(xyz0,xyz1,param)<=0.) return -1;
3b57a3f7 738 if (param[1] <= 0) {
739 param[1] = 100000000;
740 }
741 PropagateTo(x,param[1],param[0]*param[4]);
742
743 }
744
745 GetXYZ(xyz0);
746 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
747 Rotate(alpha,kTRUE);
748 GetXYZ(xyz0);
530fb4e2 749 if(GetProlongation(r,y,z)<0) return -1;
3b57a3f7 750 xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
751 xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
752 xyz1[2] = z;
753 Double_t param[7];
83dea92e 754 if(AliTracker::MeanMaterialBudget(xyz0,xyz1,param) <= 0.) return -1;
3b57a3f7 755
756 if (param[1] <= 0) {
757 param[1] = 100000000;
87a7fa94 758 }
3b57a3f7 759 PropagateTo(r,param[1],param[0]*param[4]);
760
761 return 0;
762
87a7fa94 763}
764
3b57a3f7 765//_____________________________________________________________________________
203967fc 766void AliTRDtrackV1::Print(Option_t *o) const
767{
4d6aee34 768 // Print track status
eb2b4f91 769 AliInfo(Form("PID [%4.1f %4.1f %4.1f %4.1f %4.1f]", 1.E2*fPID[0], 1.E2*fPID[1], 1.E2*fPID[2], 1.E2*fPID[3], 1.E2*fPID[4]));
203967fc 770 AliInfo(Form("Material[%5.2f %5.2f %5.2f]", fBudget[0], fBudget[1], fBudget[2]));
771
772 AliInfo(Form("x[%7.2f] t[%7.4f] alpha[%f] mass[%f]", GetX(), GetIntegratedLength(), GetAlpha(), fMass));
eb2b4f91 773 AliInfo(Form("Ntr[%1d] NtrPID[%1d] Ncl[%3d] lab[%3d]", GetNumberOfTracklets(), GetNumberOfTrackletsPID(), fN, fLab));
203967fc 774
203967fc 775 printf("|X| = (");
776 const Double_t *curP = GetParameter();
777 for (Int_t i = 0; i < 5; i++) printf("%7.2f ", curP[i]);
778 printf(")\n");
779
780 printf("|V| = \n");
781 const Double_t *curC = GetCovariance();
782 for (Int_t i = 0, j=4, k=0; i<15; i++, k++){
783 printf("%7.2f ", curC[i]);
784 if(k==j){
785 printf("\n");
786 k=-1; j--;
787 }
788 }
0217fcd0 789 if(strcmp(o, "a")!=0) return;
203967fc 790
791 for(Int_t ip=0; ip<kNplane; ip++){
792 if(!fTracklet[ip]) continue;
793 fTracklet[ip]->Print(o);
794 }
795}
796
797
798//_____________________________________________________________________________
3b57a3f7 799Bool_t AliTRDtrackV1::Rotate(Double_t alpha, Bool_t absolute)
800{
801 //
802 // Rotates track parameters in R*phi plane
803 // if absolute rotation alpha is in global system
804 // otherwise alpha rotation is relative to the current rotation angle
805 //
806
807 if (absolute) alpha -= GetAlpha();
808 //else fNRotate++;
809
810 return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
811}
87a7fa94 812
eb38ed55 813//___________________________________________________________
814void AliTRDtrackV1::SetNumberOfClusters()
815{
816// Calculate the number of clusters attached to this track
817
3b57a3f7 818 Int_t ncls = 0;
819 for(int ip=0; ip<kNplane; ip++){
68f9b6bd 820 if(fTracklet[ip] && fTrackletIndex[ip] >= 0) ncls += fTracklet[ip]->GetN();
3b57a3f7 821 }
822 AliKalmanTrack::SetNumberOfClusters(ncls);
eb38ed55 823}
824
825
0906e73e 826//_______________________________________________________________
3b57a3f7 827void AliTRDtrackV1::SetOwner()
0906e73e 828{
829 //
830 // Toggle ownership of tracklets
831 //
832
e44586fb 833 if(TestBit(kOwner)) return;
3b57a3f7 834 for (Int_t ip = 0; ip < kNplane; ip++) {
68f9b6bd 835 if(fTrackletIndex[ip]<0 || !fTracklet[ip]) continue;
3b57a3f7 836 fTracklet[ip] = new AliTRDseedV1(*fTracklet[ip]);
837 fTracklet[ip]->SetOwner();
838 }
e44586fb 839 SetBit(kOwner);
0906e73e 840}
841
842//_______________________________________________________________
4d6aee34 843void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *const trklt, Int_t index)
d9950a5a 844{
845 //
0906e73e 846 // Set the tracklets
d9950a5a 847 //
3b57a3f7 848 Int_t plane = trklt->GetPlane();
76b60503 849
3b57a3f7 850 fTracklet[plane] = trklt;
851 fTrackletIndex[plane] = index;
d9950a5a 852}
853
854//_______________________________________________________________
a310e49b 855void AliTRDtrackV1::SetTrackIn()
e3cf3d02 856{
fc6d4c2d 857// Save location of birth for the TRD track
858// If the pointer is not valid allocate memory
859//
e3cf3d02 860 const AliExternalTrackParam *op = dynamic_cast<const AliExternalTrackParam*>(this);
fc6d4c2d 861
5d4510fe 862 //printf("SetTrackIn() : fTrackLow[%p]\n", (void*)fTrackLow);
fc6d4c2d 863 if(fTrackLow){
864 fTrackLow->~AliExternalTrackParam();
865 new(fTrackLow) AliExternalTrackParam(*op);
866 } else fTrackLow = new AliExternalTrackParam(*op);
e3cf3d02 867}
868
869//_______________________________________________________________
a310e49b 870void AliTRDtrackV1::SetTrackOut(const AliExternalTrackParam *op)
e3cf3d02 871{
fc6d4c2d 872// Save location of death for the TRD track
873// If the pointer is not valid allocate memory
874//
e3cf3d02 875 if(!op) op = dynamic_cast<const AliExternalTrackParam*>(this);
fc6d4c2d 876 if(fTrackHigh){
877 fTrackHigh->~AliExternalTrackParam();
878 new(fTrackHigh) AliExternalTrackParam(*op);
879 } else fTrackHigh = new AliExternalTrackParam(*op);
e3cf3d02 880}
881
882//_______________________________________________________________
181d2c97 883void AliTRDtrackV1::UnsetTracklet(Int_t plane)
884{
87f70a90 885 if(plane<0) return;
17896e82 886 fTrackletIndex[plane] = -1;
4d6aee34 887 fTracklet[plane] = NULL;
181d2c97 888}
889
890
891//_______________________________________________________________
ef867f55 892void AliTRDtrackV1::UpdateChi2(Float_t chi2)
d9950a5a 893{
ef867f55 894// Update chi2/track with one tracklet contribution
b72f4eaf 895 SetChi2(GetChi2() + chi2);
d9950a5a 896}
897
898//_______________________________________________________________
0906e73e 899void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
d9950a5a 900{
901 //
6984f7c1 902 // Update the TRD PID information in the ESD track
d9950a5a 903 //
6984f7c1 904
93fc2389 905// Int_t nslices = AliTRDcalibDB::Instance()->GetPIDResponse(fkReconstructor->GetRecoParam()->GetPIDmethod())->GetNumberOfSlices();
e20bef2b 906 // number of tracklets used for PID calculation
907 UChar_t nPID = GetNumberOfTrackletsPID();
908 // number of tracklets attached to the track
909 UChar_t nTrk = GetNumberOfTracklets();
910 // pack the two numbers together and store them in the ESD
911 track->SetTRDntracklets(nPID | (nTrk<<3));
912 // allocate space to store raw PID signals dEdx & momentum
1f97f376 913 // independent of the method used to calculate PID (see below AliTRDPIDResponse)
914 track->SetNumberOfTRDslices((AliTRDseedV1::kNdEdxSlices+2)*AliTRDgeometry::kNlayer);
e20bef2b 915 // store raw signals
7b23a4e1 916 Float_t p, sp; Double_t spd;
6984f7c1 917 for (Int_t ip = 0; ip < kNplane; ip++) {
68f9b6bd 918 if(fTrackletIndex[ip]<0 || !fTracklet[ip]) continue;
1f97f376 919 // Fill TRD dEdx info into ESD track
920 // a. Set Summed dEdx into the first slice
921 track->SetTRDslice(fTracklet[ip]->GetdQdl(), ip, 0);
922 // b. Set NN dEdx slices
93fc2389 923 fTracklet[ip]->CookdEdx(AliTRDPIDResponse::kNslicesNN);
4d6aee34 924 const Float_t *dedx = fTracklet[ip]->GetdEdx();
1f97f376 925 for (Int_t js(0), ks(1); js < AliTRDPIDResponse::kNslicesNN; js++, ks++, dedx++){
926 if(ks>=AliTRDseedV1::kNdEdxSlices){
927 AliError(Form("Exceed allocated space for dEdx slices."));
928 break;
929 }
930 track->SetTRDslice(*dedx, ip, ks);
93fc2389 931 }
1f97f376 932 // fill TRD momentum info into ESD track
933 p = fTracklet[ip]->GetMomentum(&sp); spd = sp;
934 track->SetTRDmomentum(p, ip, &spd);
5c5d503a 935 // store global quality per tracklet instead of momentum error
936 // 26.09.11 A.Bercuci
937 // first implementation store no. of time bins filled in tracklet (5bits see "y" bits) and
938 // no. of double clusters in case of pad row cross (4bits see "x" bits)
939 // bit map for tracklet quality xxxxyyyyy
803dc399 940 // 27.10.11 A.Bercuci
c05ef2e8 941 // add chamber status bit "z" bit
803dc399 942 // bit map for tracklet quality zxxxxyyyyy
c05ef2e8 943 // 12.11.11 A.Bercuci
944 // fit tracklet quality into the field fTRDTimeBin [Char_t]
945 // bit map for tracklet quality zxxyyyyy
5554f68d 946 // The information should be retrieved by the following functions of AliESDtrack for each TRD layer
947 // GetTRDtrkltOccupancy(layer) -> no of TB filled in tracklet
948 // GetTRDtrkltClCross(layer) -> no of TB filled in crossing pad rows
949 // IsTRDtrkltChmbGood(layer) -> status of the chamber from which the tracklet is found
950 Int_t nCross(fTracklet[ip]->IsRowCross()?fTracklet[ip]->GetTBcross():0); if(nCross>3) nCross = 3;
c05ef2e8 951 Char_t trackletQ = Char_t(fTracklet[ip]->GetTBoccupancy() | (nCross<<5) | (fTracklet[ip]->IsChmbGood()<<7));
952 track->SetTRDTimBin(trackletQ, ip);
6984f7c1 953 }
e20bef2b 954 // store PID probabilities
955 track->SetTRDpid(fPID);
e2a1c98b 956
957 //store truncated mean
958 track->SetTRDsignal(fTruncatedMean);
d0037ea4 959 track->SetTRDNchamberdEdx(fNchamberdEdx);
960 track->SetTRDNclusterdEdx(fNclusterdEdx);
d9950a5a 961}
e2a1c98b 962
963//_______________________________________________________________
d0037ea4 964Double_t AliTRDtrackV1::CookTruncatedMean(const Bool_t kinvq, const Double_t mag, const Int_t charge, const Int_t kcalib, Int_t &nch, Int_t &ncls, TVectorD *Qs, TVectorD *Xs, Int_t timeBin0, Int_t timeBin1, Int_t tstep) const
e2a1c98b 965{
966 //
967 //Origin: Xianguo Lu <xianguo.lu@cern.ch>, Marian Ivanov <marian.ivanov@cern.ch>
968 //
969
970 TVectorD arrayQ(200), arrayX(200);
6951a056 971 ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, this, timeBin0, timeBin1, tstep);
e2a1c98b 972
6951a056 973 const TObjArray *cobj = kcalib ? AliTRDdEdxCalibUtils::GetObj(kinvq, mag, charge) : NULL;
e2a1c98b 974
6951a056 975 const Double_t tmean = AliTRDdEdxReconUtils::ToyCook(kinvq, ncls, &arrayQ, &arrayX, cobj);
e2a1c98b 976
6951a056 977 nch = AliTRDdEdxReconUtils::UpdateArrayX(ncls, &arrayX);
e2a1c98b 978
979 if(Qs && Xs){
980 (*Qs)=arrayQ;
981 (*Xs)=arrayX;
982 }
983
a96ab3fd 984 //printf("\ntest %.10f %d %d\n", tmean, nch, ncls);
d0037ea4 985
986 return tmean;
e2a1c98b 987}
988
24aea01c 989//_______________________________________________________________
990TObject* AliTRDtrackV1::Clone(const char* newname) const
991{
992 // temporary override TObject::Clone to avoid crashes in reco
993 AliTRDtrackV1* src = (AliTRDtrackV1*)this;
994 Bool_t isown = src->IsOwner();
f88fedf8 995 //AliInfo(Form("src_owner %d",isown));
24aea01c 996 AliTRDtrackV1* dst = new AliTRDtrackV1(*src);
997 if (isown) {
998 src->SetBit(kOwner);
999 dst->SetOwner();
1000 }
1001 return dst;
1002}