Coding rule violations corrected
[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   fNContributors(0),
40   fCovMatrix(NULL),
41   fParent(),
42   fDaughters(),
43   fProngs(NULL)
44   {
45   // default constructor
46
47   fPosition[0] = fPosition[1] = fPosition[2] = -999.;
48 }
49
50 //______________________________________________________________________________
51 AliAODVertex::AliAODVertex(const Double_t position[3], 
52                            const Double_t covMatrix[6],
53                            Double_t  chi2perNDF,
54                            TObject  *parent,
55                            Short_t id,
56                            Char_t vtype, 
57                            Int_t  nprong) :
58   AliVVertex(),
59   fChi2perNDF(chi2perNDF),
60   fID(id),
61   fType(vtype),
62   fNprong(nprong),
63   fIprong(0),
64   fNContributors(0),
65   fCovMatrix(NULL),
66   fParent(parent),
67   fDaughters(),
68   fProngs(0)
69 {
70   // constructor
71
72   SetPosition(position);
73   if (covMatrix) SetCovMatrix(covMatrix);
74   MakeProngs();
75 }
76
77 //______________________________________________________________________________
78 AliAODVertex::AliAODVertex(const Float_t position[3], 
79                            const Float_t  covMatrix[6],
80                            Double_t  chi2perNDF,
81                            TObject  *parent,
82                            Short_t id,
83                            Char_t vtype,
84                            Int_t nprong) :
85
86   AliVVertex(),
87   fChi2perNDF(chi2perNDF),
88   fID(id),
89   fType(vtype),
90   fNprong(nprong),
91   fIprong(0),
92   fNContributors(0),
93   fCovMatrix(NULL),
94   fParent(parent),
95   fDaughters(),
96   fProngs(0)
97 {
98   // constructor
99
100   SetPosition(position);
101   if (covMatrix) SetCovMatrix(covMatrix);
102   MakeProngs();
103 }
104
105 //______________________________________________________________________________
106 AliAODVertex::AliAODVertex(const Double_t position[3], 
107                            Double_t  chi2perNDF,
108                            Char_t vtype, 
109                            Int_t nprong) :
110   AliVVertex(),
111   fChi2perNDF(chi2perNDF),
112   fID(-1),
113   fType(vtype),
114   fNprong(nprong),
115   fIprong(0),
116   fNContributors(0),
117   fCovMatrix(NULL),
118   fParent(),
119   fDaughters(),
120   fProngs(0)
121 {
122   // constructor without covariance matrix
123
124   SetPosition(position);  
125   MakeProngs();
126 }
127
128 //______________________________________________________________________________
129 AliAODVertex::AliAODVertex(const Float_t position[3], 
130                            Double_t  chi2perNDF,
131                            Char_t vtype, Int_t nprong) :
132   AliVVertex(),
133   fChi2perNDF(chi2perNDF),
134   fID(-1),
135   fType(vtype),
136   fNprong(nprong),
137   fIprong(0),
138   fNContributors(0),
139   fCovMatrix(NULL),
140   fParent(),
141   fDaughters(),
142   fProngs(0)
143 {
144   // constructor without covariance matrix
145
146   SetPosition(position);  
147   MakeProngs();
148 }
149
150 //______________________________________________________________________________
151 AliAODVertex::~AliAODVertex() 
152 {
153   // Destructor
154
155   delete fCovMatrix;
156   if (fNprong > 0) delete[] fProngs;
157 }
158
159 //______________________________________________________________________________
160 AliAODVertex::AliAODVertex(const AliAODVertex& vtx) :
161   AliVVertex(vtx),
162   fChi2perNDF(vtx.fChi2perNDF),
163   fID(vtx.fID),
164   fType(vtx.fType),
165   fNprong(vtx.fNprong),
166   fIprong(vtx.fIprong),
167   fNContributors(vtx.fNContributors),
168   fCovMatrix(NULL),
169   fParent(vtx.fParent),
170   fDaughters(vtx.fDaughters),
171   fProngs(0)
172 {
173   // Copy constructor.
174   
175   for (int i = 0; i < 3; i++) 
176     fPosition[i] = vtx.fPosition[i];
177
178   if (vtx.fCovMatrix) fCovMatrix=new AliAODRedCov<3>(*vtx.fCovMatrix);
179   MakeProngs();
180   for (int i = 0; i < fNprong; i++) {
181       fProngs[i] = vtx.fProngs[i];
182   }
183 }
184
185 //______________________________________________________________________________
186 AliAODVertex* AliAODVertex::CloneWithoutRefs() const
187 {
188   // Special method to copy all but the refs 
189   
190   Double_t cov[6] = { 0.0 };
191       
192   if (fCovMatrix) fCovMatrix->GetCovMatrix(cov);
193   
194   AliAODVertex* v = new AliAODVertex(fPosition,
195                                      cov,
196                                      fChi2perNDF,
197                                      0x0,
198                                      fID,
199                                      fType,
200                                      0);
201   
202   v->SetNContributors(fNContributors);  
203   
204   return v;
205 }
206
207 //______________________________________________________________________________
208 AliAODVertex& AliAODVertex::operator=(const AliAODVertex& vtx) 
209 {
210   // Assignment operator
211   if (this != &vtx) {
212
213     // name and type
214     AliVVertex::operator=(vtx);
215
216     //momentum
217     for (int i = 0; i < 3; i++) 
218       fPosition[i] = vtx.fPosition[i];
219     
220     fChi2perNDF = vtx.fChi2perNDF;
221     fID = vtx.fID;
222     fType = vtx.fType;
223
224     //covariance matrix
225     delete fCovMatrix;
226     fCovMatrix = NULL;   
227     if (vtx.fCovMatrix) fCovMatrix=new AliAODRedCov<3>(*vtx.fCovMatrix);
228     
229     //other stuff
230     fParent = vtx.fParent;
231     fDaughters = vtx.fDaughters;
232     fNprong    = vtx.fNprong;
233     fIprong    = vtx.fIprong;  
234
235     MakeProngs();
236     for (int i = 0; i < fNprong; i++) {
237         fProngs[i] = vtx.fProngs[i];
238     }
239   }
240   
241   return *this;
242 }
243
244 //______________________________________________________________________________
245 void AliAODVertex::AddDaughter(TObject *daughter)
246 {
247   // Add reference to daughter track
248     if (!fProngs) {
249         if (fDaughters.GetEntries()==0) {
250             TRefArray* arr = &fDaughters;
251             new(arr)TRefArray(TProcessID::GetProcessWithUID(daughter));         
252         }
253         fDaughters.Add(daughter);       
254     } else {
255         if (fIprong < fNprong) {
256             fProngs[fIprong++] = daughter;
257         } else {
258             AliWarning("Number of daughters out of range !\n");
259         }
260     }
261   return;
262 }
263
264
265 //______________________________________________________________________________
266 template <class T> void AliAODVertex::GetSigmaXYZ(T sigma[3]) const
267 {
268   // Return errors on vertex position in thrust frame
269   
270   if(fCovMatrix) {
271     sigma[0]=fCovMatrix[3]; //GetCovXZ
272     sigma[1]=fCovMatrix[4]; //GetCovYZ
273     sigma[2]=fCovMatrix[5]; //GetCovZZ
274   } else 
275     sigma[0]=sigma[1]=sigma[2]=-999.;
276
277   /*
278   for (int i = 0, j = 6; i < 3; i++) {
279     j -= i+1;
280     sigma[2-i] = fCovMatrix ? TMath::Sqrt(fCovMatrix[j]) : -999.;
281   }
282   */
283 }
284
285 //______________________________________________________________________________
286 Int_t AliAODVertex::GetNContributors() const 
287 {
288   // Returns the number of tracks used to fit this vertex.
289   Int_t cont  = 0;
290
291   TString vtitle = GetTitle();
292   if (!vtitle.Contains("VertexerTracks")) {
293     cont = fNContributors;
294   } else {
295     for (Int_t iDaug = 0; iDaug < GetNDaughters(); iDaug++) {
296         AliAODTrack* aodT = dynamic_cast<AliAODTrack*>(fDaughters.At(iDaug));
297         if (!aodT) continue;
298         if (aodT->GetUsedForPrimVtxFit()) cont++;
299     } 
300     // the constraint adds another DOF
301     if(vtitle.Contains("VertexerTracksWithConstraint"))cont++;
302   }
303   return cont;
304 }
305
306 //______________________________________________________________________________
307 Bool_t AliAODVertex::HasDaughter(TObject *daughter) const 
308 {
309   // Checks if the given daughter (particle) is part of this vertex.
310     if (!fProngs) {
311         TRefArrayIter iter(&fDaughters);
312         while (TObject *daugh = iter.Next()) {
313             if (daugh == daughter) return kTRUE;
314         }
315         return kFALSE;
316     } else {
317         Bool_t has = kFALSE;
318         for (int i = 0; i < fNprong; i++) {
319             if (fProngs[i].GetObject() == daughter) has = kTRUE;
320         }
321         return has;
322     }
323 }
324
325 //______________________________________________________________________________
326 Double_t AliAODVertex::RotatedCovMatrixXX(Double_t phi, Double_t theta) const
327 {
328   // XX term of covariance matrix after rotation by phi around z-axis
329   // and, then, by theta around new y-axis
330
331   if (!fCovMatrix) {
332     //AliFatal("Covariance matrix not set");
333     return -999.;
334   }
335
336   Double_t covMatrix[6];
337
338   GetCovMatrix(covMatrix);
339
340   Double_t cp = TMath::Cos(phi);
341   Double_t sp = TMath::Sin(phi);
342   Double_t ct = TMath::Cos(theta);
343   Double_t st = TMath::Sin(theta);
344   return
345      covMatrix[0]*cp*cp*ct*ct  // GetCovXX
346     +covMatrix[1]*2.*cp*sp*ct*ct  // GetCovXY
347     +covMatrix[3]*2.*cp*ct*st  // GetCovXZ
348     +covMatrix[2]*sp*sp*ct*ct  // GetCovYY
349     +covMatrix[4]*2.*sp*ct*st  // GetCovYZ
350     +covMatrix[5]*st*st;  // GetCovZZ
351 }
352
353 //______________________________________________________________________________
354 Double_t AliAODVertex::RotatedCovMatrixXY(Double_t phi, Double_t theta) const
355 {
356   // XY term of covariance matrix after rotation by phi around z-axis
357   // and, then, by theta around new y-axis
358
359   if (!fCovMatrix) {
360     //AliFatal("Covariance matrix not set");
361     return -999.;
362   }
363
364   Double_t covMatrix[6];
365
366   GetCovMatrix(covMatrix);
367
368   Double_t cp = TMath::Cos(phi);
369   Double_t sp = TMath::Sin(phi);
370   Double_t ct = TMath::Cos(theta);
371   Double_t st = TMath::Sin(theta);
372   return 
373     -covMatrix[0]*cp*sp*ct  // GetCovXX
374     +covMatrix[1]*ct*(cp*cp-sp*sp)  // GetCovXY
375     -covMatrix[3]*sp*st  // GetCovXZ
376     +covMatrix[2]*cp*sp*ct  // GetCovYY
377     +covMatrix[4]*cp*st;  // GetCovYZ
378 }
379
380 //______________________________________________________________________________
381 Double_t AliAODVertex::RotatedCovMatrixXZ(Double_t phi, Double_t theta) const
382 {
383   // XZ term of covariance matrix after rotation by phi around z-axis
384   // and, then, by theta around new y-axis
385
386   if (!fCovMatrix) {
387     //AliFatal("Covariance matrix not set");
388     return -999.;
389   }
390
391   Double_t covMatrix[6];
392
393   GetCovMatrix(covMatrix);
394
395   Double_t cp = TMath::Cos(phi);
396   Double_t sp = TMath::Sin(phi);
397   Double_t ct = TMath::Cos(theta);
398   Double_t st = TMath::Sin(theta);
399   return 
400     -covMatrix[0]*cp*cp*ct*st  // GetCovXX
401     -covMatrix[1]*2.*cp*sp*ct*st  // GetCovXY
402     +covMatrix[3]*cp*(ct*ct-st*st)  // GetCovXZ
403     -covMatrix[2]*sp*sp*ct*st  // GetCovYY
404     +covMatrix[4]*sp*(ct*ct-st*st)  // GetCovYZ
405     +covMatrix[5]*ct*st;  // GetCovZZ
406 }
407
408 //______________________________________________________________________________
409 Double_t AliAODVertex::RotatedCovMatrixYY(Double_t phi) const
410 {
411   // YY term of covariance matrix after rotation by phi around z-axis
412   // and, then, by theta around new y-axis
413
414   if (!fCovMatrix) {
415     //AliFatal("Covariance matrix not set");
416     return -999.;
417   }
418
419   Double_t covMatrix[6];
420
421   GetCovMatrix(covMatrix);
422
423   Double_t cp = TMath::Cos(phi);
424   Double_t sp = TMath::Sin(phi);
425   return
426      covMatrix[0]*sp*sp  // GetCovXX
427     -covMatrix[1]*2.*cp*sp  // GetCovXY
428     +covMatrix[2]*cp*cp;  // GetCovYY
429 }
430
431 //______________________________________________________________________________
432 Double_t AliAODVertex::RotatedCovMatrixYZ(Double_t phi, Double_t theta) const
433 {
434   // YZ term of covariance matrix after rotation by phi around z-axis
435   // and, then, by theta around new y-axis
436
437   if (!fCovMatrix) {
438     //AliFatal("Covariance matrix not set");
439     return -999.;
440   }
441
442   Double_t covMatrix[6];
443
444   GetCovMatrix(covMatrix);
445
446   Double_t cp = TMath::Cos(phi);
447   Double_t sp = TMath::Sin(phi);
448   Double_t ct = TMath::Cos(theta);
449   Double_t st = TMath::Sin(theta);
450   return 
451      covMatrix[0]*cp*sp*st  // GetCovXX
452     +covMatrix[1]*st*(sp*sp-cp*cp)  // GetCovXY
453     -covMatrix[3]*sp*ct  // GetCovXZ
454     -covMatrix[2]*cp*sp*st  // GetCovYY
455     +covMatrix[4]*cp*ct;  // GetCovYZ
456 }
457
458 //______________________________________________________________________________
459 Double_t AliAODVertex::RotatedCovMatrixZZ(Double_t phi, Double_t theta) const
460 {
461   // ZZ term of covariance matrix after rotation by phi around z-axis
462   // and, then, by theta around new y-axis
463
464   if (!fCovMatrix) {
465     //AliFatal("Covariance matrix not set");
466     return -999.;
467   }
468
469   Double_t covMatrix[6];
470
471   GetCovMatrix(covMatrix);
472
473   Double_t cp = TMath::Cos(phi);
474   Double_t sp = TMath::Sin(phi);
475   Double_t ct = TMath::Cos(theta);
476   Double_t st = TMath::Sin(theta);
477   return
478      covMatrix[0]*cp*cp*st*st  // GetCovXX
479     +covMatrix[1]*2.*cp*sp*st*st  // GetCovXY
480     -covMatrix[3]*2.*cp*ct*st  // GetCovXZ
481     +covMatrix[2]*sp*sp*st*st  // GetCovYY
482     -covMatrix[4]*2.*sp*sp*ct*st  // GetCovYZ
483     +covMatrix[5]*ct*ct;  // GetCovZZ
484 }
485
486 //______________________________________________________________________________
487 Double_t AliAODVertex::Distance2ToVertex(const AliAODVertex *vtx) const
488 {
489   // distance in 3D to another AliAODVertex
490
491   Double_t dx = GetX()-vtx->GetX();
492   Double_t dy = GetY()-vtx->GetY();
493   Double_t dz = GetZ()-vtx->GetZ();
494
495   return dx*dx+dy*dy+dz*dz;
496 }
497
498 //______________________________________________________________________________
499 Double_t AliAODVertex::DistanceXY2ToVertex(const AliAODVertex *vtx) const
500 {
501   // distance in XY to another AliAODVertex
502
503   Double_t dx = GetX()-vtx->GetX();
504   Double_t dy = GetY()-vtx->GetY();
505
506   return dx*dx+dy*dy;
507 }
508
509 //______________________________________________________________________________
510 Double_t AliAODVertex::Error2DistanceToVertex(AliAODVertex *vtx) const
511 {
512   // error on the distance in 3D to another AliAODVertex
513
514   Double_t phi,theta;
515   PhiAndThetaToVertex(vtx,phi,theta);
516   // error2 due to this vertex
517   Double_t error2 = RotatedCovMatrixXX(phi,theta);
518   // error2 due to vtx vertex
519   Double_t error2vtx = vtx->RotatedCovMatrixXX(phi,theta);
520
521   return error2+error2vtx;
522 }
523
524 //______________________________________________________________________________
525 Double_t AliAODVertex::Error2DistanceXYToVertex(AliAODVertex *vtx) const
526 {
527   // error on the distance in XY to another AliAODVertex
528
529   Double_t phi,theta;
530   PhiAndThetaToVertex(vtx,phi,theta);
531   // error2 due to this vertex
532   Double_t error2 = RotatedCovMatrixXX(phi);
533   // error2 due to vtx vertex
534   Double_t error2vtx = vtx->RotatedCovMatrixXX(phi);
535
536   return error2+error2vtx;
537 }
538
539 //______________________________________________________________________________
540 template <class T, class P> 
541 void AliAODVertex::PhiAndThetaToVertex(AliAODVertex *vtx, P &phi, T &theta) const
542 {
543   // rotation angles around z-axis (phi) and around new y-axis (theta)
544   // with which vtx is seen (used by RotatedCovMatrix... methods)
545
546   phi = TMath::Pi()+TMath::ATan2(-vtx->GetY()+GetY(),-vtx->GetX()+GetX());
547   Double_t vtxxphi = vtx->GetX()*TMath::Cos(phi)+vtx->GetY()*TMath::Sin(phi);
548   Double_t xphi = GetX()*TMath::Cos(phi)+GetY()*TMath::Sin(phi);
549   theta = TMath::ATan2(vtx->GetZ()-GetZ(),vtxxphi-xphi);
550 }
551
552 //______________________________________________________________________________
553 void AliAODVertex::PrintIndices() const 
554 {
555   // Print indices of particles originating form this vertex
556
557   TRefArrayIter iter(&fDaughters);
558   while (TObject *daugh = iter.Next()) {
559     printf("Particle %p originates from this vertex.\n", static_cast<void*>(daugh));
560   }
561 }
562
563 //______________________________________________________________________________
564 const char* AliAODVertex::AsString() const
565 {
566   // Make a string describing this object
567   
568   TString tmp(Form("%10s pos(%7.2f,%7.2f,%7.2f)",GetTypeName((AODVtx_t)GetType()),GetX(),GetY(),GetZ()));
569   
570   if (GetType()==kPrimary || GetType()==kMainSPD || GetType()==kPileupSPD )
571   {
572     tmp += Form(" ncontrib %d chi2/ndf %4.1f",GetNContributors(),GetChi2perNDF());
573
574   }
575   
576   if ( !fParent.GetObject() ) 
577   {
578     tmp += " no parent";
579   }
580   if ( fDaughters.GetEntriesFast() > 0 )
581   {
582     if ( fDaughters.GetEntriesFast() == 1 ) 
583     {
584       tmp += " origin of 1 particle";
585     }
586     else
587     {
588       tmp += Form(" origin of %2d particles",fDaughters.GetEntriesFast());
589     }
590   }
591   
592   return tmp.Data();
593 }
594
595 //______________________________________________________________________________
596 const char* AliAODVertex::GetTypeName(AODVtx_t type)
597 {
598   // Return an ASCII version of type
599   
600   switch (type)
601   {
602     case kPrimary:
603       return "primary";
604       break;
605     case kKink:
606       return "kink";
607       break;
608     case kV0:
609       return "v0";
610       break;
611     case kCascade:
612       return "cascade";
613       break;
614     case kMainSPD:
615       return "mainSPD";
616       break;
617     case kPileupSPD:
618       return "pileupSPD";
619       break;
620     case kPileupTracks:
621       return "pileupTRK";
622       break;
623     case kMainTPC:
624       return "mainTPC";
625       break;
626     default:
627       return "unknown";
628       break;
629   };
630 }
631
632 //______________________________________________________________________________
633 void AliAODVertex::Print(Option_t* /*option*/) const 
634 {
635   // Print information of all data members
636
637   printf("Vertex position:\n");
638   printf("     x = %f\n", fPosition[0]);
639   printf("     y = %f\n", fPosition[1]);
640   printf("     z = %f\n", fPosition[2]);
641   printf(" parent particle: %p\n", static_cast<void*>(fParent.GetObject()));
642   printf(" origin of %d particles\n", fDaughters.GetEntriesFast());
643   printf(" vertex type %d\n", fType);
644   
645   /*
646   if (fCovMatrix) {
647     printf("Covariance matrix:\n");
648     printf(" %12.10f  %12.10f  %12.10f\n %12.10f  %12.10f  %12.10f\n %12.10f  %12.10f  %12.10f\n", 
649            fCovMatrix[0],
650            fCovMatrix[1],
651            fCovMatrix[3],
652            fCovMatrix[1],
653            fCovMatrix[2],
654            fCovMatrix[4],
655            fCovMatrix[3],
656            fCovMatrix[4],
657            fCovMatrix[5]); 
658            } */
659   printf(" Chi^2/NDF = %f\n", fChi2perNDF);
660 }
661