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