remove obselete clean QA macros command
[u/mrichter/AliRoot.git] / PWG1 / 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 }
85
86
87 AliESDRecInfo::AliESDRecInfo(const AliESDRecInfo& recinfo):
88   TObject(),
89   fITSOn(0),           // ITS refitted inward
90   fTRDOn(0),           // ITS refitted inward
91   fDeltaP(0),          //delta of momenta
92   fSign(0),           // sign
93   fReconstructed(0),         //flag if track was reconstructed
94   fFake(0),             // fake track
95   fMultiple(0),         // number of reconstructions
96   fTPCOn(0),           // TPC refitted inward
97   fBestTOFmatch(0),        //best matching between times
98   fESDtrack(0),        // esd track
99   fTrackF(0),      // friend track
100   fTPCtrack(0),        // tpc track
101   fITStrack(0),        // its track
102   fTRDtrack(0),        // trd track  
103   fTracks(0)           // array of tracks with the same label
104 {
105   //
106   //
107   //
108   memcpy(this,&recinfo, sizeof(recinfo));
109   fESDtrack=0; fTrackF=0; fTPCtrack=0;fITStrack=0;fTRDtrack=0;
110   fTracks=0;
111   SetESDtrack(recinfo.GetESDtrack());
112 }
113
114
115 AliESDRecInfo& AliESDRecInfo::operator=(const AliESDRecInfo& info) { 
116   //
117   // Assignment operator
118   //
119   this->~AliESDRecInfo();
120   new (this) AliESDRecInfo(info);
121   return *this;
122 }
123
124
125
126 AliESDRecInfo::~AliESDRecInfo()
127
128 {
129   //
130   //  destructor
131   //
132   if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
133   if (fTrackF)   {   fTrackF=0;}
134   if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
135   if (fITStrack) { delete fITStrack; fITStrack=0;}
136   if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
137   if (fTracks) { 
138     delete  fTracks;  fTracks=0;
139   }
140 }
141
142
143
144 void AliESDRecInfo::Reset()
145 {
146   //
147   // reset info
148   //
149   fMultiple =0; 
150   fFake     =0;
151   fReconstructed=0;
152   if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
153   if (fTrackF)   { fTrackF=0;}
154   if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
155   if (fITStrack) { delete fITStrack; fITStrack=0;}
156   if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
157   if (fTracks) { delete  fTracks;  fTracks=0;}
158
159
160 void AliESDRecInfo::SetESDtrack(const AliESDtrack *track){
161   //
162   // 
163   //
164
165   if (fESDtrack) delete fESDtrack;
166   fESDtrack = (AliESDtrack*)track->Clone();
167   //AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(fESDtrack->GetID());
168
169
170   if (track->GetFriendTrack()){
171     fTrackF = (AliESDfriendTrack*)track->GetFriendTrack();
172     Int_t icalib=0;
173     TObject *cobject=0;
174     //
175     while (fTrackF->GetCalibObject(icalib)){
176       cobject=fTrackF->GetCalibObject(icalib);
177       if (dynamic_cast<AliTPCseed*>(cobject)){
178         if (fTPCtrack) delete fTPCtrack;
179         fTPCtrack = (AliTPCseed*)(dynamic_cast<AliTPCseed*>(cobject))->Clone();
180       } else if (dynamic_cast<AliTRDtrackV1*>(cobject)){
181         if (fTRDtrack) delete fTRDtrack;
182         fTRDtrack = (AliTRDtrackV1*)(dynamic_cast<AliTRDtrackV1*>(cobject))->Clone();
183       }
184       icalib++;
185     }
186   }
187   if (!fTPCtrack) fTPCtrack = new AliTPCseed; // add dummy track
188 }
189
190
191
192 void AliESDRecInfo::AddESDtrack(const AliESDtrack *track, AliMCInfo* info){
193   //
194   // Add ESDtrack
195   //
196   AliESDtrack *nctrack = (AliESDtrack*) track;
197   fMultiple++;
198   if (!fESDtrack) {
199     SetESDtrack(track);  
200     Update(info,0,kTRUE);
201     return;
202   }
203   if (!fTracks) fTracks = new TClonesArray("AliESDtrack",10);
204   Int_t ntracks = fTracks->GetEntriesFast();
205   new ((*fTracks)[ntracks]) AliESDtrack(*track);
206   if  (nctrack->GetKinkIndex(0)>0) return; 
207   //
208   //
209   //
210   if (!nctrack->IsOn(AliESDtrack::kTPCrefit) && fStatus[1]==2) return;
211   if (!nctrack->IsOn(AliESDtrack::kITSin) && fStatus[0]>0) return;
212   if ( nctrack->GetParameter()[4]*info->GetCharge()<0) return; //returning track 
213   
214   Float_t dtheta = TMath::ATan(nctrack->GetTgl())-info->GetParticle().Theta()-TMath::Pi()-2;
215   if (TMath::Abs(dtheta)>0.1) return;
216
217   SetESDtrack(track);  
218   Update(info,0,kTRUE);
219 }
220
221
222
223 void  AliESDRecInfo::UpdatePoints(AliESDtrack*track)
224 {
225   //
226   //
227   Int_t iclusters[200];
228   Float_t density[160];
229   for (Int_t i=0;i<160;i++) density[i]=-1.;
230   fTPCPoints[0]= 160;
231   fTPCPoints[1] = -1;
232   //
233   if (fTPCPoints[0]<fTPCPoints[1]) return;
234   //  Int_t nclusters=track->GetTPCclusters(iclusters);
235
236   Int_t ngood=0;
237   Int_t undeff=0;
238   Int_t nall =0;
239   Int_t range=20;
240   for (Int_t i=0;i<160;i++){
241     Int_t last = i-range;
242     if (nall<range) nall++;
243     if (last>=0){
244       if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--;
245       if (iclusters[last]==-1) undeff--;
246     }
247     if (iclusters[i]>0&& (iclusters[i]&0x8000)==0)   ngood++;
248     if (iclusters[i]==-1) undeff++;
249     if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
250   }
251   Float_t maxdens=0;
252   Int_t indexmax =0;
253   for (Int_t i=0;i<160;i++){
254     if (density[i]<0) continue;
255     if (density[i]>maxdens){
256       maxdens=density[i];
257       indexmax=i;
258     }
259   }
260   //
261   //max dens point
262   fTPCPoints[3] = maxdens;
263   fTPCPoints[1] = indexmax;
264   //
265   // last point
266   for (Int_t i=indexmax;i<160;i++){
267     if (density[i]<0) continue;
268     if (density[i]<maxdens/2.) {
269       break;
270     }
271     fTPCPoints[2]=i;
272   }
273   //
274   // first point
275   for (Int_t i=indexmax;i>0;i--){
276     if (density[i]<0) continue;
277     if (density[i]<maxdens/2.) {
278       break;
279     }
280     fTPCPoints[0]=i;
281   }
282   //
283   // Density at the last 30 padrows
284   //
285   // 
286   nall  = 0;
287   ngood = 0;
288   for (Int_t i=159;i>0;i--){
289     if (iclusters[i]==-1) continue; //dead zone
290     nall++;
291     if (iclusters[i]>0)   ngood++;
292     if (nall>20) break;
293   }
294   fTPCPoints[4] = Float_t(ngood)/Float_t(nall);
295   //
296   if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1;
297
298
299 }
300
301 //
302 //
303 void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed)
304 {
305   //
306   //calculates derived variables
307   //  
308   UpdatePoints(fESDtrack);
309   UpdateStatus(info,reconstructed);
310   UpdateITS(info);
311   UpdateTPC(info);
312   UpdateTOF(info);
313 }
314
315
316 void  AliESDRecInfo::UpdateStatus(AliMCInfo* info, Bool_t reconstructed){
317   //
318   // Interpret bit mask flags
319   //
320   for (Int_t i=0;i<4;i++) fStatus[i] =0;
321   fReconstructed = kFALSE;
322   fTPCOn = kFALSE;
323   fITSOn = kFALSE;
324   fTRDOn = kFALSE;  
325   if (reconstructed==kFALSE) return;
326
327   fLabels[0] = info->fLabel;
328   fLabels[1] = info->fPrimPart;
329   fReconstructed = kTRUE;
330   fTPCOn = ((fESDtrack->GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
331   fITSOn = ((fESDtrack->GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
332   fTRDOn = ((fESDtrack->GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
333   //
334   //  
335   if ((fESDtrack->GetStatus()&AliESDtrack::kTPCrefit)>0){
336     fStatus[1] =3;
337   }
338   else{
339     if ((fESDtrack->GetStatus()&AliESDtrack::kTPCout)>0){
340       fStatus[1] =2;
341     }
342     else{
343       if ((fESDtrack->GetStatus()&AliESDtrack::kTPCin)>0)
344         fStatus[1]=1;
345     }      
346   }
347   //
348   if ((fESDtrack->GetStatus()&AliESDtrack::kITSout)>0){
349     fStatus[0] =2;
350   }
351   else{
352     if ((fESDtrack->GetStatus()&AliESDtrack::kITSrefit)>0){
353       fStatus[0] =1;
354     }
355     else{
356       fStatus[0]=0;
357     }      
358   }
359   //
360   //
361   if ((fESDtrack->GetStatus()&AliESDtrack::kTRDrefit)>0){
362     fStatus[2] =2;
363   }
364   else{
365     if ((fESDtrack->GetStatus()&AliESDtrack::kTRDout)>0){
366       fStatus[2] =1;
367     }
368   }
369   if ((fESDtrack->GetStatus()&AliESDtrack::kTRDStop)>0){
370     fStatus[2] =10;
371   }   
372 }
373
374 void AliESDRecInfo::UpdateTOF(AliMCInfo* info){
375   //
376   // Update TOF related comparison information
377   //
378   fBestTOFmatch=1000;
379   if (((fESDtrack->GetStatus()&AliESDtrack::kTOFout)>0)){
380     //
381     // best tof match
382     Double_t times[5];
383     fESDtrack->GetIntegratedTimes(times);    
384     for (Int_t i=0;i<5;i++){
385       if ( TMath::Abs(fESDtrack->GetTOFsignal()-times[i]) <TMath::Abs(fBestTOFmatch) ){
386         fBestTOFmatch = fESDtrack->GetTOFsignal()-times[i];
387       }
388     }
389     Int_t toflabel[3];
390     fESDtrack->GetTOFLabel(toflabel);
391     Bool_t toffake=kTRUE;
392     Bool_t tofdaughter=kFALSE;
393     for (Int_t i=0;i<3;i++){
394       if (toflabel[i]<0) continue;      
395       if (toflabel[i]== TMath::Abs(fESDtrack->GetLabel()))  toffake=kFALSE;     
396       if (toflabel[i]==info->fParticle.GetDaughter(0) || (toflabel[i]==info->fParticle.GetDaughter(1))) tofdaughter=kTRUE;  // decay product of original particle
397       fStatus[3]=1;
398     }
399     if (toffake) fStatus[3] =3;       //total fake
400     if (tofdaughter) fStatus[3]=2;    //fake because of decay
401   }else{
402     fStatus[3]=0;
403   }
404 }
405
406
407
408 void AliESDRecInfo::UpdateITS(AliMCInfo* info){
409   //
410   // Update ITS related comparison information
411   //
412   fITSinP0[0]=info->fParticle.Px();
413   fITSinP0[1]=info->fParticle.Py();
414   fITSinP0[2]=info->fParticle.Pz();
415   fITSinP0[3]=info->fParticle.Pt();    
416   //
417   fITSinR0[0]=info->fParticle.Vx();
418   fITSinR0[1]=info->fParticle.Vy();
419   fITSinR0[2]=info->fParticle.Vz();
420   fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
421   fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
422   //
423   //
424   if (fITSinP0[3]>0.0000001){
425     fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
426     fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
427   }
428   if (fITSOn){
429     // ITS 
430     Double_t param[5],x;
431     fESDtrack->GetExternalParameters(x,param);   
432     //    fESDtrack->GetConstrainedExternalParameters(x,param);   
433     Double_t cov[15];
434     fESDtrack->GetExternalCovariance(cov);
435     //fESDtrack->GetConstrainedExternalCovariance(cov);
436     if (TMath::Abs(param[4])<0.0000000001) return;
437     
438     fESDtrack->GetXYZ(fITSinR1);
439     fESDtrack->GetPxPyPz(fITSinP1);
440     fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]);
441     //
442     fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]);
443     fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]);
444     //
445     //
446     if (fITSinP1[3]>0.0000001){
447       fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]);
448       fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]);  
449     }
450     //
451     //
452     fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3];  //delta rfi
453     fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]);
454     fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]);              //delta z
455     fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]);
456     fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]);
457     fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]);
458     fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1]));
459     fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
460     fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);    
461     Double_t sign = (param[4]>0) ? 1:-1; 
462     fSign = sign;
463     fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);
464   }
465 }
466
467 void AliESDRecInfo::UpdateTPC(AliMCInfo* info){
468   //
469   //
470   // Update TPC related information
471   //
472   AliTrackReference * ref = &(info->fTrackRef);
473   fTPCinR0[0] = info->fTrackRef.X();    
474   fTPCinR0[1] = info->fTrackRef.Y();    
475   fTPCinR0[2] = info->fTrackRef.Z();
476   fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
477   fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
478   //
479   fTPCinP0[0] = ref->Px();
480   fTPCinP0[1] = ref->Py();
481   fTPCinP0[2] = ref->Pz();
482   fTPCinP0[3] = ref->Pt();
483   fTPCinP0[4] = ref->P();
484   fDeltaP     = (ref->P()-info->fParticle.P())/info->fParticle.P();
485   //
486   //
487   if (fTPCinP0[3]>0.0000001){
488     //
489     fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
490     fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
491   }
492   //
493   if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){
494     //TPC
495     fESDtrack->GetInnerXYZ(fTPCinR1);
496     fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
497     fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);        
498     fESDtrack->GetInnerPxPyPz(fTPCinP1);
499     fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
500     fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]);
501     //
502     //
503     if (fTPCinP1[3]>0.000000000000001){
504       fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
505       fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);  
506     }    
507     Double_t cov[15], param[5],x, alpha;
508     fESDtrack->GetInnerExternalCovariance(cov);
509     fESDtrack->GetInnerExternalParameters(alpha, x,param);
510     if (x<50) return ;
511     //
512     fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3];  //delta rfi
513     fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
514     fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]);              //delta z
515     fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
516     fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
517     fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
518     fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]));
519     fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
520     fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
521     Double_t sign = (param[4]>0)? 1.:-1; 
522     fSign =sign;
523     fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(TMath::Abs(cov[14]));
524   }
525 }
526
527
528
529
530