9ee1d1f8f1bc06237849dfc776624349513c5cec
[u/mrichter/AliRoot.git] / STEER / AliAODVertex.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, 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 /* $Id$ */
17
18 //-------------------------------------------------------------------------
19 //     AOD track base class
20 //     Base class for Analysis Object Data
21 //     Generic version
22 //     Author: Markus Oldenburg, CERN
23 //     Inheritance from AliVVertex: A. Dainese
24 //-------------------------------------------------------------------------
25
26 #include "AliAODVertex.h"
27 #include "AliAODTrack.h"
28
29 ClassImp(AliAODVertex)
30
31 //______________________________________________________________________________
32 AliAODVertex::AliAODVertex() : 
33   AliVVertex(),
34   fChi2perNDF(-999.),
35   fID(-1),
36   fType(kUndef),
37   fNprong(0),
38   fIprong(0),
39   fCovMatrix(NULL),
40   fParent(),
41   fDaughters(),
42   fProngs(NULL)
43   {
44   // default constructor
45
46   fPosition[0] = fPosition[1] = fPosition[2] = -999.;
47 }
48
49 //______________________________________________________________________________
50 AliAODVertex::AliAODVertex(const Double_t position[3], 
51                            const Double_t covMatrix[6],
52                            Double_t  chi2perNDF,
53                            TObject  *parent,
54                            Short_t id,
55                            Char_t vtype, 
56                            Int_t  nprong) :
57   AliVVertex(),
58   fChi2perNDF(chi2perNDF),
59   fID(id),
60   fType(vtype),
61   fNprong(nprong),
62   fIprong(0),
63   fCovMatrix(NULL),
64   fParent(parent),
65   fDaughters(),
66   fProngs(0)
67 {
68   // constructor
69
70   SetPosition(position);
71   if (covMatrix) SetCovMatrix(covMatrix);
72   MakeProngs();
73 }
74
75 //______________________________________________________________________________
76 AliAODVertex::AliAODVertex(const Float_t position[3], 
77                            const Float_t  covMatrix[6],
78                            Double_t  chi2perNDF,
79                            TObject  *parent,
80                            Short_t id,
81                            Char_t vtype,
82                            Int_t nprong) :
83
84   AliVVertex(),
85   fChi2perNDF(chi2perNDF),
86   fID(id),
87   fType(vtype),
88   fNprong(nprong),
89   fIprong(0),
90   fCovMatrix(NULL),
91   fParent(parent),
92   fDaughters(),
93   fProngs(0)
94 {
95   // constructor
96
97   SetPosition(position);
98   if (covMatrix) SetCovMatrix(covMatrix);
99   MakeProngs();
100 }
101
102 //______________________________________________________________________________
103 AliAODVertex::AliAODVertex(const Double_t position[3], 
104                            Double_t  chi2perNDF,
105                            Char_t vtype, 
106                            Int_t nprong) :
107   AliVVertex(),
108   fChi2perNDF(chi2perNDF),
109   fID(-1),
110   fType(vtype),
111   fNprong(nprong),
112   fIprong(0),
113   fCovMatrix(NULL),
114   fParent(),
115   fDaughters(),
116   fProngs(0)
117 {
118   // constructor without covariance matrix
119
120   SetPosition(position);  
121   MakeProngs();
122 }
123
124 //______________________________________________________________________________
125 AliAODVertex::AliAODVertex(const Float_t position[3], 
126                            Double_t  chi2perNDF,
127                            Char_t vtype, Int_t nprong) :
128   AliVVertex(),
129   fChi2perNDF(chi2perNDF),
130   fID(-1),
131   fType(vtype),
132   fNprong(nprong),
133   fIprong(0),
134   fCovMatrix(NULL),
135   fParent(),
136   fDaughters(),
137   fProngs(0)
138 {
139   // constructor without covariance matrix
140
141   SetPosition(position);  
142   MakeProngs();
143 }
144
145 //______________________________________________________________________________
146 AliAODVertex::~AliAODVertex() 
147 {
148   // Destructor
149
150   delete fCovMatrix;
151   if (fNprong > 0) delete[] fProngs;
152 }
153
154 //______________________________________________________________________________
155 AliAODVertex::AliAODVertex(const AliAODVertex& vtx) :
156   AliVVertex(vtx),
157   fChi2perNDF(vtx.fChi2perNDF),
158   fID(vtx.fID),
159   fType(vtx.fType),
160   fNprong(vtx.fNprong),
161   fIprong(vtx.fIprong),
162   fCovMatrix(NULL),
163   fParent(vtx.fParent),
164   fDaughters(vtx.fDaughters),
165   fProngs(0)
166 {
167   // Copy constructor.
168   
169   for (int i = 0; i < 3; i++) 
170     fPosition[i] = vtx.fPosition[i];
171
172   if (vtx.fCovMatrix) fCovMatrix=new AliAODRedCov<3>(*vtx.fCovMatrix);
173   MakeProngs();
174   for (int i = 0; i < fNprong; i++) {
175       fProngs[i] = vtx.fProngs[i];
176   }
177 }
178
179 //______________________________________________________________________________
180 AliAODVertex& AliAODVertex::operator=(const AliAODVertex& vtx) 
181 {
182   // Assignment operator
183   if (this != &vtx) {
184
185     // name and type
186     AliVVertex::operator=(vtx);
187
188     //momentum
189     for (int i = 0; i < 3; i++) 
190       fPosition[i] = vtx.fPosition[i];
191     
192     fChi2perNDF = vtx.fChi2perNDF;
193     fID = vtx.fID;
194     fType = vtx.fType;
195
196     //covariance matrix
197     delete fCovMatrix;
198     fCovMatrix = NULL;   
199     if (vtx.fCovMatrix) fCovMatrix=new AliAODRedCov<3>(*vtx.fCovMatrix);
200     
201     //other stuff
202     fParent = vtx.fParent;
203     fDaughters = vtx.fDaughters;
204     fNprong    = vtx.fNprong;
205     fIprong    = vtx.fIprong;  
206
207     MakeProngs();
208     for (int i = 0; i < fNprong; i++) {
209         fProngs[i] = vtx.fProngs[i];
210     }
211   }
212   
213   return *this;
214 }
215
216 //______________________________________________________________________________
217 void AliAODVertex::AddDaughter(TObject *daughter)
218 {
219   // Add reference to daughter track
220     if (!fProngs) {
221         if (fDaughters.GetEntries()==0) {
222             TRefArray* arr = &fDaughters;
223             new(arr)TRefArray(TProcessID::GetProcessWithUID(daughter));         
224         }
225         fDaughters.Add(daughter);       
226     } else {
227         if (fIprong < fNprong) {
228             fProngs[fIprong++] = daughter;
229         } else {
230             AliWarning("Number of daughters out of range !\n");
231         }
232     }
233   return;
234 }
235
236
237 //______________________________________________________________________________
238 template <class T> void AliAODVertex::GetSigmaXYZ(T sigma[3]) const
239 {
240   // Return errors on vertex position in thrust frame
241   
242   if(fCovMatrix) {
243     sigma[0]=fCovMatrix[3]; //GetCovXZ
244     sigma[1]=fCovMatrix[4]; //GetCovYZ
245     sigma[2]=fCovMatrix[5]; //GetCovZZ
246   } else 
247     sigma[0]=sigma[1]=sigma[2]=-999.;
248
249   /*
250   for (int i = 0, j = 6; i < 3; i++) {
251     j -= i+1;
252     sigma[2-i] = fCovMatrix ? TMath::Sqrt(fCovMatrix[j]) : -999.;
253   }
254   */
255 }
256
257 //______________________________________________________________________________
258 Int_t AliAODVertex::GetNContributors() const 
259 {
260   // Returns the number of tracks used to fit this vertex.
261   Int_t cont = 0;
262
263   if (!strcmp(GetTitle(), "vertexer: 3D")) {
264     cont = fNContributors;
265   } else {
266     for (Int_t iDaug = 0; iDaug < GetNDaughters(); iDaug++) {
267       if (((AliAODTrack*)fDaughters.At(iDaug))->GetUsedForVtxFit()) cont++;
268     }
269   }
270   return cont;
271 }
272
273 //______________________________________________________________________________
274 Bool_t AliAODVertex::HasDaughter(TObject *daughter) const 
275 {
276   // Checks if the given daughter (particle) is part of this vertex.
277     if (!fProngs) {
278         TRefArrayIter iter(&fDaughters);
279         while (TObject *daugh = iter.Next()) {
280             if (daugh == daughter) return kTRUE;
281         }
282         return kFALSE;
283     } else {
284         Bool_t has = kFALSE;
285         for (int i = 0; i < fNprong; i++) {
286             if (fProngs[i].GetObject() == daughter) has = kTRUE;
287         }
288         return has;
289     }
290 }
291
292 //______________________________________________________________________________
293 Double_t AliAODVertex::RotatedCovMatrixXX(Double_t phi, Double_t theta) const
294 {
295   // XX term of covariance matrix after rotation by phi around z-axis
296   // and, then, by theta around new y-axis
297
298   if (!fCovMatrix) {
299     //AliFatal("Covariance matrix not set");
300     return -999.;
301   }
302
303   Double_t covMatrix[6];
304
305   GetCovMatrix(covMatrix);
306
307   Double_t cp = TMath::Cos(phi);
308   Double_t sp = TMath::Sin(phi);
309   Double_t ct = TMath::Cos(theta);
310   Double_t st = TMath::Sin(theta);
311   return
312      covMatrix[0]*cp*cp*ct*ct  // GetCovXX
313     +covMatrix[1]*2.*cp*sp*ct*ct  // GetCovXY
314     +covMatrix[3]*2.*cp*ct*st  // GetCovXZ
315     +covMatrix[2]*sp*sp*ct*ct  // GetCovYY
316     +covMatrix[4]*2.*sp*ct*st  // GetCovYZ
317     +covMatrix[5]*st*st;  // GetCovZZ
318 }
319
320 //______________________________________________________________________________
321 Double_t AliAODVertex::RotatedCovMatrixXY(Double_t phi, Double_t theta) const
322 {
323   // XY term of covariance matrix after rotation by phi around z-axis
324   // and, then, by theta around new y-axis
325
326   if (!fCovMatrix) {
327     //AliFatal("Covariance matrix not set");
328     return -999.;
329   }
330
331   Double_t covMatrix[6];
332
333   GetCovMatrix(covMatrix);
334
335   Double_t cp = TMath::Cos(phi);
336   Double_t sp = TMath::Sin(phi);
337   Double_t ct = TMath::Cos(theta);
338   Double_t st = TMath::Sin(theta);
339   return 
340     -covMatrix[0]*cp*sp*ct  // GetCovXX
341     +covMatrix[1]*ct*(cp*cp-sp*sp)  // GetCovXY
342     -covMatrix[3]*sp*st  // GetCovXZ
343     +covMatrix[2]*cp*sp*ct  // GetCovYY
344     +covMatrix[4]*cp*st;  // GetCovYZ
345 }
346
347 //______________________________________________________________________________
348 Double_t AliAODVertex::RotatedCovMatrixXZ(Double_t phi, Double_t theta) const
349 {
350   // XZ term of covariance matrix after rotation by phi around z-axis
351   // and, then, by theta around new y-axis
352
353   if (!fCovMatrix) {
354     //AliFatal("Covariance matrix not set");
355     return -999.;
356   }
357
358   Double_t covMatrix[6];
359
360   GetCovMatrix(covMatrix);
361
362   Double_t cp = TMath::Cos(phi);
363   Double_t sp = TMath::Sin(phi);
364   Double_t ct = TMath::Cos(theta);
365   Double_t st = TMath::Sin(theta);
366   return 
367     -covMatrix[0]*cp*cp*ct*st  // GetCovXX
368     -covMatrix[1]*2.*cp*sp*ct*st  // GetCovXY
369     +covMatrix[3]*cp*(ct*ct-st*st)  // GetCovXZ
370     -covMatrix[2]*sp*sp*ct*st  // GetCovYY
371     +covMatrix[4]*sp*(ct*ct-st*st)  // GetCovYZ
372     +covMatrix[5]*ct*st;  // GetCovZZ
373 }
374
375 //______________________________________________________________________________
376 Double_t AliAODVertex::RotatedCovMatrixYY(Double_t phi) const
377 {
378   // YY term of covariance matrix after rotation by phi around z-axis
379   // and, then, by theta around new y-axis
380
381   if (!fCovMatrix) {
382     //AliFatal("Covariance matrix not set");
383     return -999.;
384   }
385
386   Double_t covMatrix[6];
387
388   GetCovMatrix(covMatrix);
389
390   Double_t cp = TMath::Cos(phi);
391   Double_t sp = TMath::Sin(phi);
392   return
393      covMatrix[0]*sp*sp  // GetCovXX
394     -covMatrix[1]*2.*cp*sp  // GetCovXY
395     +covMatrix[2]*cp*cp;  // GetCovYY
396 }
397
398 //______________________________________________________________________________
399 Double_t AliAODVertex::RotatedCovMatrixYZ(Double_t phi, Double_t theta) const
400 {
401   // YZ term of covariance matrix after rotation by phi around z-axis
402   // and, then, by theta around new y-axis
403
404   if (!fCovMatrix) {
405     //AliFatal("Covariance matrix not set");
406     return -999.;
407   }
408
409   Double_t covMatrix[6];
410
411   GetCovMatrix(covMatrix);
412
413   Double_t cp = TMath::Cos(phi);
414   Double_t sp = TMath::Sin(phi);
415   Double_t ct = TMath::Cos(theta);
416   Double_t st = TMath::Sin(theta);
417   return 
418      covMatrix[0]*cp*sp*st  // GetCovXX
419     +covMatrix[1]*st*(sp*sp-cp*cp)  // GetCovXY
420     -covMatrix[3]*sp*ct  // GetCovXZ
421     -covMatrix[2]*cp*sp*st  // GetCovYY
422     +covMatrix[4]*cp*ct;  // GetCovYZ
423 }
424
425 //______________________________________________________________________________
426 Double_t AliAODVertex::RotatedCovMatrixZZ(Double_t phi, Double_t theta) const
427 {
428   // ZZ term of covariance matrix after rotation by phi around z-axis
429   // and, then, by theta around new y-axis
430
431   if (!fCovMatrix) {
432     //AliFatal("Covariance matrix not set");
433     return -999.;
434   }
435
436   Double_t covMatrix[6];
437
438   GetCovMatrix(covMatrix);
439
440   Double_t cp = TMath::Cos(phi);
441   Double_t sp = TMath::Sin(phi);
442   Double_t ct = TMath::Cos(theta);
443   Double_t st = TMath::Sin(theta);
444   return
445      covMatrix[0]*cp*cp*st*st  // GetCovXX
446     +covMatrix[1]*2.*cp*sp*st*st  // GetCovXY
447     -covMatrix[3]*2.*cp*ct*st  // GetCovXZ
448     +covMatrix[2]*sp*sp*st*st  // GetCovYY
449     -covMatrix[4]*2.*sp*sp*ct*st  // GetCovYZ
450     +covMatrix[5]*ct*ct;  // GetCovZZ
451 }
452
453 //______________________________________________________________________________
454 Double_t AliAODVertex::DistanceToVertex(AliAODVertex *vtx) const
455 {
456   // distance in 3D to another AliAODVertex
457
458   Double_t dx = GetX()-vtx->GetX();
459   Double_t dy = GetY()-vtx->GetY();
460   Double_t dz = GetZ()-vtx->GetZ();
461
462   return TMath::Sqrt(dx*dx+dy*dy+dz*dz);
463 }
464
465 //______________________________________________________________________________
466 Double_t AliAODVertex::DistanceXYToVertex(AliAODVertex *vtx) const
467 {
468   // distance in XY to another AliAODVertex
469
470   Double_t dx = GetX()-vtx->GetX();
471   Double_t dy = GetY()-vtx->GetY();
472
473   return TMath::Sqrt(dx*dx+dy*dy);
474 }
475
476 //______________________________________________________________________________
477 Double_t AliAODVertex::ErrorDistanceToVertex(AliAODVertex *vtx) const
478 {
479   // error on the distance in 3D to another AliAODVertex
480
481   Double_t phi,theta;
482   PhiAndThetaToVertex(vtx,phi,theta);
483   // error2 due to this vertex
484   Double_t error2 = RotatedCovMatrixXX(phi,theta);
485   // error2 due to vtx vertex
486   Double_t error2vtx = vtx->RotatedCovMatrixXX(phi,theta);
487
488   return TMath::Sqrt(error2+error2vtx);
489 }
490
491 //______________________________________________________________________________
492 Double_t AliAODVertex::ErrorDistanceXYToVertex(AliAODVertex *vtx) const
493 {
494   // error on the distance in XY to another AliAODVertex
495
496   Double_t phi,theta;
497   PhiAndThetaToVertex(vtx,phi,theta);
498   // error2 due to this vertex
499   Double_t error2 = RotatedCovMatrixXX(phi);
500   // error2 due to vtx vertex
501   Double_t error2vtx = vtx->RotatedCovMatrixXX(phi);
502
503   return TMath::Sqrt(error2+error2vtx);
504 }
505
506 //______________________________________________________________________________
507 template <class T, class P> 
508 void AliAODVertex::PhiAndThetaToVertex(AliAODVertex *vtx, P &phi, T &theta) const
509 {
510   // rotation angles around z-axis (phi) and around new y-axis (theta)
511   // with which vtx is seen (used by RotatedCovMatrix... methods)
512
513   phi = TMath::Pi()+TMath::ATan2(-vtx->GetY()+GetY(),-vtx->GetX()+GetX());
514   Double_t vtxxphi = vtx->GetX()*TMath::Cos(phi)+vtx->GetY()*TMath::Sin(phi);
515   Double_t xphi = GetX()*TMath::Cos(phi)+GetY()*TMath::Sin(phi);
516   theta = TMath::ATan2(vtx->GetZ()-GetZ(),vtxxphi-xphi);
517 }
518
519 //______________________________________________________________________________
520 void AliAODVertex::PrintIndices() const 
521 {
522   // Print indices of particles originating form this vertex
523
524   TRefArrayIter iter(&fDaughters);
525   while (TObject *daugh = iter.Next()) {
526     printf("Particle %p originates from this vertex.\n", static_cast<void*>(daugh));
527   }
528 }
529
530 //______________________________________________________________________________
531 void AliAODVertex::Print(Option_t* /*option*/) const 
532 {
533   // Print information of all data members
534
535   printf("Vertex position:\n");
536   printf("     x = %f\n", fPosition[0]);
537   printf("     y = %f\n", fPosition[1]);
538   printf("     z = %f\n", fPosition[2]);
539   printf(" parent particle: %p\n", static_cast<void*>(fParent.GetObject()));
540   printf(" origin of %d particles\n", fDaughters.GetEntriesFast());
541   printf(" vertex type %d\n", fType);
542   
543   /*
544   if (fCovMatrix) {
545     printf("Covariance matrix:\n");
546     printf(" %12.10f  %12.10f  %12.10f\n %12.10f  %12.10f  %12.10f\n %12.10f  %12.10f  %12.10f\n", 
547            fCovMatrix[0],
548            fCovMatrix[1],
549            fCovMatrix[3],
550            fCovMatrix[1],
551            fCovMatrix[2],
552            fCovMatrix[4],
553            fCovMatrix[3],
554            fCovMatrix[4],
555            fCovMatrix[5]); 
556            } */
557   printf(" Chi^2/NDF = %f\n", fChi2perNDF);
558 }
559