STEER/AOD/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