don't lie in the log!
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliESDRecInfo.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17///////////////////////////////////////////////////////////////////////////////
18// //
19// Time Projection Chamber //
20// Comparison macro for reconstructed tracks //
21// responsible:
22// marian.ivanov@cern.ch //
23//
24//
25
26
27
28
29
30//ROOT includes
31#include "Rtypes.h"
32//
33//ALIROOT includes
34//
35#include "AliESDtrack.h"
36#include "AliTracker.h"
37#include "AliTPCParam.h"
38#include "AliTrackReference.h"
39#include "AliTPCParamSR.h"
40#include "AliESDfriend.h"
41#include "AliESDtrack.h"
42#include "AliTPCseed.h"
43#include "AliITStrackMI.h"
44#include "AliTRDtrackV1.h"
45#include "AliMCInfo.h"
46#include "AliESDRecInfo.h"
47
48
49
50ClassImp(AliESDRecInfo)
51
52
53
54
55AliTPCParam * GetTPCParam(){
56 AliTPCParamSR * par = new AliTPCParamSR;
57 par->Update();
58 return par;
59}
60
61
62
63
64AliESDRecInfo::AliESDRecInfo():
65 fITSOn(0), // ITS refitted inward
66 fTRDOn(0), // ITS refitted inward
67 fDeltaP(0), //delta of momenta
68 fSign(0), // sign
69 fReconstructed(0), //flag if track was reconstructed
70 fFake(0), // fake track
71 fMultiple(0), // number of reconstructions
72 fTPCOn(0), // TPC refitted inward
73 fBestTOFmatch(0), //best matching between times
74 fESDtrack(0), // esd track
75 fTrackF(0), // friend track
76 fTPCtrack(0), // tpc track
77 fITStrack(0), // its track
78 fTRDtrack(0), // trd track
79 fTracks(0) // array of tracks with the same label
80{
81 //
82 // default constructor
83 //
84 for(Int_t i=0; i<10; i++) { fTPCPoints[i] = 0;}
85 for(Int_t i=0; i<5; i++) {
86 fTPCinR0[i] = 0 ;
87 fTPCinR1[i] = 0 ;
88 fTPCinP0[i] = 0 ;
89 fTPCinP1[i] = 0 ;
90 fTPCDelta[i] = 0 ;
91 fTPCPools[i] = 0 ;
92 fITSinR0[i] = 0 ;
93 fITSinR1[i] = 0 ;
94 fITSinP0[i] = 0 ;
95 fITSinP1[i] = 0 ;
96 fITSDelta[i] = 0 ;
97 fITSPools[i] = 0 ;
98 }
99 for(Int_t i=0; i<2; i++) {
100 fTPCAngle0[i] = 0 ;
101 fTPCAngle1[i] = 0 ;
102 fITSAngle0[i] = 0 ;
103 fITSAngle1[i] = 0 ;
104 fLabels[i] = 0 ;
105 }
106 for(Int_t i=0; i<3; i++) { fTRLocalCoord[i] = 0 ; }
107 for(Int_t i=0; i<4; i++) { fStatus[i] = 0 ; }
108
109}
110
111
112AliESDRecInfo::AliESDRecInfo(const AliESDRecInfo& recinfo):
113 TObject(),
114 fITSOn(0), // ITS refitted inward
115 fTRDOn(0), // ITS refitted inward
116 fDeltaP(0), //delta of momenta
117 fSign(0), // sign
118 fReconstructed(0), //flag if track was reconstructed
119 fFake(0), // fake track
120 fMultiple(0), // number of reconstructions
121 fTPCOn(0), // TPC refitted inward
122 fBestTOFmatch(0), //best matching between times
123 fESDtrack(0), // esd track
124 fTrackF(0), // friend track
125 fTPCtrack(0), // tpc track
126 fITStrack(0), // its track
127 fTRDtrack(0), // trd track
128 fTracks(0) // array of tracks with the same label
129{
130 //
131 //
132 //
133 for(Int_t i=0; i<10; i++) { fTPCPoints[i] = 0;}
134 for(Int_t i=0; i<5; i++) {
135 fTPCinR0[i] = 0 ;
136 fTPCinR1[i] = 0 ;
137 fTPCinP0[i] = 0 ;
138 fTPCinP1[i] = 0 ;
139 fTPCDelta[i] = 0 ;
140 fTPCPools[i] = 0 ;
141 fITSinR0[i] = 0 ;
142 fITSinR1[i] = 0 ;
143 fITSinP0[i] = 0 ;
144 fITSinP1[i] = 0 ;
145 fITSDelta[i] = 0 ;
146 fITSPools[i] = 0 ;
147 }
148 for(Int_t i=0; i<2; i++) {
149 fTPCAngle0[i] = 0 ;
150 fTPCAngle1[i] = 0 ;
151 fITSAngle0[i] = 0 ;
152 fITSAngle1[i] = 0 ;
153 fLabels[i] = 0 ;
154 }
155 for(Int_t i=0; i<3; i++) { fTRLocalCoord[i] = 0 ; }
156 for(Int_t i=0; i<4; i++) { fStatus[i] = 0 ; }
157
158 memcpy(this,&recinfo, sizeof(recinfo));
159 fESDtrack=0; fTrackF=0; fTPCtrack=0;fITStrack=0;fTRDtrack=0;
160 fTracks=0;
161 SetESDtrack(recinfo.GetESDtrack());
162
163}
164
165
166AliESDRecInfo& AliESDRecInfo::operator=(const AliESDRecInfo& info) {
167 //
168 // Assignment operator
169 //
170 this->~AliESDRecInfo();
171 new (this) AliESDRecInfo(info);
172 return *this;
173}
174
175
176
177AliESDRecInfo::~AliESDRecInfo()
178
179{
180 //
181 // destructor
182 //
183 if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
184 if (fTrackF) { fTrackF=0;}
185 if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
186 if (fITStrack) { delete fITStrack; fITStrack=0;}
187 if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
188 if (fTracks) {
189 delete fTracks; fTracks=0;
190 }
191}
192
193
194
195void AliESDRecInfo::Reset()
196{
197 //
198 // reset info
199 //
200 fMultiple =0;
201 fFake =0;
202 fReconstructed=0;
203 if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
204 if (fTrackF) { fTrackF=0;}
205 if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
206 if (fITStrack) { delete fITStrack; fITStrack=0;}
207 if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
208 if (fTracks) { delete fTracks; fTracks=0;}
209}
210
211void AliESDRecInfo::SetESDtrack(const AliESDtrack *track){
212 //
213 //
214 //
215
216 if (fESDtrack) delete fESDtrack;
217 fESDtrack = (AliESDtrack*)track->Clone();
218 //AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(fESDtrack->GetID());
219
220
221 if (track->GetFriendTrack()){
222 fTrackF = (AliESDfriendTrack*)track->GetFriendTrack();
223 Int_t icalib=0;
224 TObject *cobject=0;
225 //
226 AliTPCseed *tpcSeed = NULL;
227 AliTRDtrackV1 *trdTrack = NULL;
228 while (fTrackF->GetCalibObject(icalib)){
229 cobject=fTrackF->GetCalibObject(icalib);
230 if ((tpcSeed = dynamic_cast<AliTPCseed*>(cobject))){
231 if (fTPCtrack) delete fTPCtrack;
232 fTPCtrack = (AliTPCseed*)tpcSeed->Clone();
233 } else if ((trdTrack = dynamic_cast<AliTRDtrackV1*>(cobject))){
234 if (fTRDtrack) delete fTRDtrack;
235 fTRDtrack = (AliTRDtrackV1*)trdTrack->Clone();
236 }
237 icalib++;
238 }
239 }
240 if (!fTPCtrack) fTPCtrack = new AliTPCseed; // add dummy track
241}
242
243
244
245void AliESDRecInfo::AddESDtrack(const AliESDtrack *track, AliMCInfo* info){
246 //
247 // Add ESDtrack
248 //
249 AliESDtrack *nctrack = (AliESDtrack*) track;
250 fMultiple++;
251 if (!fESDtrack) {
252 SetESDtrack(track);
253 Update(info,0,kTRUE);
254 return;
255 }
256 if (!fTracks) fTracks = new TClonesArray("AliESDtrack",10);
257 Int_t ntracks = fTracks->GetEntriesFast();
258 new ((*fTracks)[ntracks]) AliESDtrack(*track);
259 if (nctrack->GetKinkIndex(0)>0) return;
260 //
261 //
262 //
263 if (!nctrack->IsOn(AliESDtrack::kTPCrefit) && fStatus[1]==2) return;
264 if (!nctrack->IsOn(AliESDtrack::kITSin) && fStatus[0]>0) return;
265 if ( nctrack->GetParameter()[4]*info->GetCharge()<0) return; //returning track
266
267 Float_t dtheta = TMath::ATan(nctrack->GetTgl())-info->GetParticle().Theta()-TMath::Pi()-2;
268 if (TMath::Abs(dtheta)>0.1) return;
269
270 SetESDtrack(track);
271 Update(info,0,kTRUE);
272}
273
274
275
276void AliESDRecInfo::UpdatePoints(AliESDtrack*track)
277{
278 //
279 //
280 Int_t iclusters[200];
281 Float_t density[160];
282 for (Int_t i=0;i<160;i++) density[i]=-1.;
283 fTPCPoints[0]= 160;
284 fTPCPoints[1] = -1;
285 //
286 if (fTPCPoints[0]<fTPCPoints[1]) return;
287 /*Int_t nclusters=*/track->GetTPCclusters(iclusters);
288
289 Int_t ngood=0;
290 Int_t undeff=0;
291 Int_t nall =0;
292 Int_t range=20;
293 for (Int_t i=0;i<160;i++){
294 Int_t last = i-range;
295 if (nall<range) nall++;
296 if (last>=0){
297 if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--;
298 if (iclusters[last]==-1) undeff--;
299 }
300 if (iclusters[i]>0&& (iclusters[i]&0x8000)==0) ngood++;
301 if (iclusters[i]==-1) undeff++;
302 if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
303 }
304 Float_t maxdens=0;
305 Int_t indexmax =0;
306 for (Int_t i=0;i<160;i++){
307 if (density[i]<0) continue;
308 if (density[i]>maxdens){
309 maxdens=density[i];
310 indexmax=i;
311 }
312 }
313 //
314 //max dens point
315 fTPCPoints[3] = maxdens;
316 fTPCPoints[1] = indexmax;
317 //
318 // last point
319 for (Int_t i=indexmax;i<160;i++){
320 if (density[i]<0) continue;
321 if (density[i]<maxdens/2.) {
322 break;
323 }
324 fTPCPoints[2]=i;
325 }
326 //
327 // first point
328 for (Int_t i=indexmax;i>0;i--){
329 if (density[i]<0) continue;
330 if (density[i]<maxdens/2.) {
331 break;
332 }
333 fTPCPoints[0]=i;
334 }
335 //
336 // Density at the last 30 padrows
337 //
338 //
339 nall = 0;
340 ngood = 0;
341 for (Int_t i=159;i>0;i--){
342 if (iclusters[i]==-1) continue; //dead zone
343 nall++;
344 if (iclusters[i]>0) ngood++;
345 if (nall>20) break;
346 }
347 fTPCPoints[4] = Float_t(ngood)/Float_t(nall);
348 //
349 if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1;
350
351
352}
353
354//
355//
356void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed)
357{
358 //
359 //calculates derived variables
360 //
361 UpdatePoints(fESDtrack);
362 UpdateStatus(info,reconstructed);
363 UpdateITS(info);
364 UpdateTPC(info);
365 UpdateTOF(info);
366}
367
368
369void AliESDRecInfo::UpdateStatus(AliMCInfo* info, Bool_t reconstructed){
370 //
371 // Interpret bit mask flags
372 //
373 for (Int_t i=0;i<4;i++) fStatus[i] =0;
374 fReconstructed = kFALSE;
375 fTPCOn = kFALSE;
376 fITSOn = kFALSE;
377 fTRDOn = kFALSE;
378 if (reconstructed==kFALSE) return;
379
380 fLabels[0] = info->fLabel;
381 fLabels[1] = info->fPrimPart;
382 fReconstructed = kTRUE;
383 fTPCOn = ((fESDtrack->GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
384 fITSOn = ((fESDtrack->GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
385 fTRDOn = ((fESDtrack->GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
386 //
387 //
388 if ((fESDtrack->GetStatus()&AliESDtrack::kTPCrefit)>0){
389 fStatus[1] =3;
390 }
391 else{
392 if ((fESDtrack->GetStatus()&AliESDtrack::kTPCout)>0){
393 fStatus[1] =2;
394 }
395 else{
396 if ((fESDtrack->GetStatus()&AliESDtrack::kTPCin)>0)
397 fStatus[1]=1;
398 }
399 }
400 //
401 if ((fESDtrack->GetStatus()&AliESDtrack::kITSout)>0){
402 fStatus[0] =2;
403 }
404 else{
405 if ((fESDtrack->GetStatus()&AliESDtrack::kITSrefit)>0){
406 fStatus[0] =1;
407 }
408 else{
409 fStatus[0]=0;
410 }
411 }
412 //
413 //
414 if ((fESDtrack->GetStatus()&AliESDtrack::kTRDrefit)>0){
415 fStatus[2] =2;
416 }
417 else{
418 if ((fESDtrack->GetStatus()&AliESDtrack::kTRDout)>0){
419 fStatus[2] =1;
420 }
421 }
422 if ((fESDtrack->GetStatus()&AliESDtrack::kTRDStop)>0){
423 fStatus[2] =10;
424 }
425}
426
427void AliESDRecInfo::UpdateTOF(AliMCInfo* info){
428 //
429 // Update TOF related comparison information
430 //
431 fBestTOFmatch=1000;
432 if (((fESDtrack->GetStatus()&AliESDtrack::kTOFout)>0)){
433 //
434 // best tof match
435 Double_t times[AliPID::kSPECIESC];
436 fESDtrack->GetIntegratedTimes(times,AliPID::kSPECIESC);
437 for (Int_t i=0;i<AliPID::kSPECIESC;i++){
438 if ( TMath::Abs(fESDtrack->GetTOFsignal()-times[i]) <TMath::Abs(fBestTOFmatch) ){
439 fBestTOFmatch = fESDtrack->GetTOFsignal()-times[i];
440 }
441 }
442 Int_t toflabel[3];
443 fESDtrack->GetTOFLabel(toflabel);
444 Bool_t toffake=kTRUE;
445 Bool_t tofdaughter=kFALSE;
446 for (Int_t i=0;i<3;i++){
447 if (toflabel[i]<0) continue;
448 if (toflabel[i]== TMath::Abs(fESDtrack->GetLabel())) toffake=kFALSE;
449 if (toflabel[i]==info->fParticle.GetDaughter(0) || (toflabel[i]==info->fParticle.GetDaughter(1))) tofdaughter=kTRUE; // decay product of original particle
450 fStatus[3]=1;
451 }
452 if (toffake) fStatus[3] =3; //total fake
453 if (tofdaughter) fStatus[3]=2; //fake because of decay
454 }else{
455 fStatus[3]=0;
456 }
457}
458
459
460
461void AliESDRecInfo::UpdateITS(AliMCInfo* info){
462 //
463 // Update ITS related comparison information
464 //
465 fITSinP0[0]=info->fParticle.Px();
466 fITSinP0[1]=info->fParticle.Py();
467 fITSinP0[2]=info->fParticle.Pz();
468 fITSinP0[3]=info->fParticle.Pt();
469 //
470 fITSinR0[0]=info->fParticle.Vx();
471 fITSinR0[1]=info->fParticle.Vy();
472 fITSinR0[2]=info->fParticle.Vz();
473 fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
474 fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
475 //
476 //
477 if (fITSinP0[3]>0.0000001){
478 fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
479 fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
480 }
481 if (fITSOn){
482 // ITS
483 Double_t param[5],x;
484 fESDtrack->GetExternalParameters(x,param);
485 // fESDtrack->GetConstrainedExternalParameters(x,param);
486 Double_t cov[15];
487 fESDtrack->GetExternalCovariance(cov);
488 //fESDtrack->GetConstrainedExternalCovariance(cov);
489 if (TMath::Abs(param[4])<0.0000000001) return;
490
491 fESDtrack->GetXYZ(fITSinR1);
492 fESDtrack->GetPxPyPz(fITSinP1);
493 fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]);
494 //
495 fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]);
496 fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]);
497 //
498 //
499 if (fITSinP1[3]>0.0000001){
500 fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]);
501 fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]);
502 }
503 //
504 //
505 fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3]; //delta rfi
506 fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]);
507 fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]); //delta z
508 fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]);
509 fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]);
510 fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]);
511 fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1]));
512 fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
513 fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);
514 Double_t sign = (param[4]>0) ? 1:-1;
515 fSign = sign;
516 fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);
517 }
518}
519
520void AliESDRecInfo::UpdateTPC(AliMCInfo* info){
521 //
522 //
523 // Update TPC related information
524 //
525 AliTrackReference * ref = &(info->fTrackRef);
526 fTPCinR0[0] = info->fTrackRef.X();
527 fTPCinR0[1] = info->fTrackRef.Y();
528 fTPCinR0[2] = info->fTrackRef.Z();
529 fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
530 fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
531 //
532 fTPCinP0[0] = ref->Px();
533 fTPCinP0[1] = ref->Py();
534 fTPCinP0[2] = ref->Pz();
535 fTPCinP0[3] = ref->Pt();
536 fTPCinP0[4] = ref->P();
537 fDeltaP = (ref->P()-info->fParticle.P())/info->fParticle.P();
538 //
539 //
540 if (fTPCinP0[3]>0.0000001){
541 //
542 fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
543 fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
544 }
545 //
546 if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){
547 //TPC
548 fESDtrack->GetInnerXYZ(fTPCinR1);
549 fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
550 fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);
551 fESDtrack->GetInnerPxPyPz(fTPCinP1);
552 fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
553 fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]);
554 //
555 //
556 if (fTPCinP1[3]>0.000000000000001){
557 fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
558 fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);
559 }
560 Double_t cov[15], param[5],x, alpha;
561 fESDtrack->GetInnerExternalCovariance(cov);
562 fESDtrack->GetInnerExternalParameters(alpha, x,param);
563 if (x<50) return ;
564 //
565 fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3]; //delta rfi
566 fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
567 fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]); //delta z
568 fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
569 fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
570 fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
571 fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]));
572 fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
573 fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
574 Double_t sign = (param[4]>0)? 1.:-1;
575 fSign =sign;
576 fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(TMath::Abs(cov[14]));
577 }
578}
579
580
581
582
583