don't lie in the log!
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliESDRecInfo.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 ///////////////////////////////////////////////////////////////////////////////
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
50 ClassImp(AliESDRecInfo)
51
52
53
54
55 AliTPCParam * GetTPCParam(){
56   AliTPCParamSR * par = new AliTPCParamSR;
57   par->Update();
58   return par;
59 }
60
61
62
63
64 AliESDRecInfo::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
112 AliESDRecInfo::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
166 AliESDRecInfo& 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
177 AliESDRecInfo::~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
195 void 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
211 void 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
245 void 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
276 void  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 //
356 void 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
369 void  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
427 void 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
461 void 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
520 void 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