RALICE/AliTrack.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 // $Id$
  17
  18 ///////////////////////////////////////////////////////////////////////////
  19 // Class AliTrack
  20 // Handling of the attributes of a reconstructed particle track.
  21 //
  22 // Coding example :
  23 // ----------------
  24 //
  25 // Float_t a[4]={195.,1.2,-0.04,8.5};
  26 // Ali4Vector pmu;
  27 // pmu.SetVector(a,"car");
  28 // AliTrack t1;
  29 // t1.Set4Momentum(pmu);
  30 //
  31 // Float_t b[3]={1.2,-0.04,8.5};
  32 // Ali3Vector p;
  33 // p.SetVector(b,"car");
  34 // AliTrack t2;
  35 // t2.Set3Momentum(p);
  36 // t2.SetCharge(0);
  37 // t2.SetMass(1.115);
  38 //
  39 // t1.Info();
  40 // t2.Info();
  41 //
  42 // Float_t pi=acos(-1.);
  43 // Float_t thcms=0.2*pi; // decay theta angle in cms
  44 // Float_t phicms=pi/4.; // decay theta angle in cms
  45 // Float_t m1=0.938;
  46 // Float_t m2=0.140;
  47 // t2.Decay(m1,m2,thcms,phicms); // Track t2 decay : Lambda -> proton + pion
  48 //
  49 // t2.List();
  50 //
  51 // Int_t ndec=t2.GetNdecay();
  52 // AliTrack* d1=t2.GetDecayTrack(1); // Access to decay track number 1
  53 // AliTrack* d2=t2.GetDecayTrack(2); // Access to decay track number 2
  54 //
  55 // AliSignal s1,s2,s3,s4;
  56 //
  57 // .... // Code (e.g. detector readout) to fill AliSignal data
  58 //
  59 // AliTrack trec; // Track which will be reconstructed from signals
  60 // trec.AddSignal(s1);
  61 // trec.AddSignal(s3);
  62 // trec.AddSignal(s4);
  63 //
  64 // Ali3Vector P;
  65 // Float_t Q,M;
  66 //
  67 // ... // Code which accesses signals from trec and reconstructs
  68 //        3-momentum P, charge Q, mass M etc...
  69 //
  70 // trec.Set3Momentum(P);
  71 // trec.SetCharge(Q);
  72 // trec.SetMass(M);
  73 //
  74 // Float_t r1[3]={1.6,-3.8,25.7};
  75 // Float_t er1[3]={0.2,0.5,1.8};
  76 // Float_t r2[3]={8.6,23.8,-6.7};
  77 // Float_t er2[3]={0.93,1.78,0.8};
  78 // AliPosition begin,end;
  79 // begin.SetPosition(r1,"car");
  80 // begin.SetPositionErrors(er1,"car");
  81 // end.SetPosition(r2,"car");
  82 // end.SetPositionErrors(er2,"car");
  83 // trec.SetBeginPoint(begin);
  84 // trec.SetEndPoint(end);
  85 //
  86 // Note : All quantities are in GeV, GeV/c or GeV/c**2
  87 //
  88 //--- Author: Nick van Eijndhoven 10-jul-1997 UU-SAP Utrecht
  89 //- Modified: NvE $Date$ UU-SAP Utrecht
  90 ///////////////////////////////////////////////////////////////////////////
  91
  92 #include "AliTrack.h"
  93
  94 ClassImp(AliTrack) // Class implementation to enable ROOT I/O
  95
  96 AliTrack::AliTrack()
  97 {
  98 // Default constructor
  99 // All variables initialised to 0
 100  fDecays=0;
 101  fSignals=0;
 102  fMasses=0;
 103  fDmasses=0;
 104  fPmasses=0;
 105  Reset();
 106 }
 107 ///////////////////////////////////////////////////////////////////////////
 108 AliTrack::~AliTrack()
 109 {
 110 // Destructor to delete memory allocated for decay tracks array
 111  if (fDecays)
 112  {
 113   fDecays->Delete();
 114   delete fDecays;
 115   fDecays=0;
 116  }
 117  if (fSignals)
 118  {
 119   fSignals->Clear();
 120   delete fSignals;
 121   fSignals=0;
 122  }
 123 }
 124 ///////////////////////////////////////////////////////////////////////////
 125 void AliTrack::Reset()
 126 {
 127 // Reset all variables to 0
 128  fQ=0;
 129  fNdec=0;
 130  fNsig=0;
 131  fNmasses=0;
 132  Double_t a[4]={0,0,0,0};
 133  SetVector(a,"sph");
 134  if (fDecays)
 135  {
 136   fDecays->Delete();
 137   delete fDecays;
 138   fDecays=0;
 139  }
 140  if (fSignals)
 141  {
 142   fSignals->Clear();
 143   delete fSignals;
 144   fSignals=0;
 145  }
 146  Double_t b[3]={0,0,0};
 147  fBegin.SetPosition(b,"sph");
 148  fEnd.SetPosition(b,"sph");
 149  if (fMasses)
 150  {
 151   delete fMasses;
 152   fMasses=0;
 153  }
 154  if (fDmasses)
 155  {
 156   delete fDmasses;
 157   fDmasses=0;
 158  }
 159  if (fPmasses)
 160  {
 161   delete fPmasses;
 162   fPmasses=0;
 163  }
 164 }
 165 ///////////////////////////////////////////////////////////////////////////
 166 void AliTrack::Set3Momentum(Ali3Vector& p)
 167 {
 168 // Set the track parameters according to the 3-momentum p
 169  Set3Vector(p);
 170 }
 171 ///////////////////////////////////////////////////////////////////////////
 172 void AliTrack::Set4Momentum(Ali4Vector& p)
 173 {
 174 // Set the track parameters according to the 4-momentum p
 175  Double_t E=p.GetScalar();
 176  Double_t dE=p.GetResultError();
 177  Ali3Vector pv=p.Get3Vector();
 178  SetVector(E,pv);
 179  SetScalarError(dE);
 180 }
 181 ///////////////////////////////////////////////////////////////////////////
 182 void AliTrack::SetMass(Double_t m,Double_t dm)
 183 {
 184 // Set the particle mass
 185 // The default value for the error dm is 0.
 186  Double_t inv=pow(m,2);
 187  Double_t dinv=fabs(2.*m*dm);
 188  SetInvariant(inv,dinv);
 189 }
 190 ///////////////////////////////////////////////////////////////////////////
 191 void AliTrack::SetCharge(Float_t q)
 192 {
 193 // Set the particle charge
 194  fQ=q;
 195 }
 196 ///////////////////////////////////////////////////////////////////////////
 197 void AliTrack::Info(TString f)
 198 {
 199 // Provide track information within the coordinate frame f
 200  Double_t m=GetMass();
 201  Double_t dm=GetResultError();
 202  cout << " *AliTrack::Info* Mass : " << m
 203       << " error : " << dm << " Charge : " << fQ
 204       << " Momentum : " << GetMomentum() << " Nmass hyp. : " << fNmasses
 205       << " Ntracks : " << fNdec << " Nsignals : " << fNsig << endl;
 206  for (Int_t i=0; i<fNmasses; i++)
 207  {
 208   cout << " Mass hypothesis " << (i+1) << " Mass : " << fMasses->At(i)
 209        << " error : " << fDmasses->At(i) << " prob. : " << fPmasses->At(i)
 210        << endl;
 211  }
 212  Ali4Vector::Info(f);
 213 }
 214 ///////////////////////////////////////////////////////////////////////////
 215 void AliTrack::List(TString f)
 216 {
 217 // Provide current track and decay level 1 information within coordinate frame f
 218
 219  Info(f); // Information of the current track
 220
 221  // Decay products of this track
 222  AliTrack* td;
 223  for (Int_t id=1; id<=fNdec; id++)
 224  {
 225   td=GetDecayTrack(id);
 226   if (td)
 227   {
 228    cout << "  ---Level 1 sec. track no. " << id << endl;
 229    td->Info(f);
 230   }
 231   else
 232   {
 233    cout << " *AliTrack::List* Error : No decay track present." << endl;
 234   }
 235  }
 236 }
 237 ///////////////////////////////////////////////////////////////////////////
 238 void AliTrack::ListAll(TString f)
 239 {
 240 // Provide complete track and decay information within the coordinate frame f
 241
 242  Info(f); // Information of the current track
 243  cout << " Begin-point :"; fBegin.Info(f);
 244  cout << " End-point   :"; fEnd.Info(f);
 245  for (Int_t is=1; is<=GetNsignals(); is++)
 246  {
 247   ((AliSignal*)GetSignal(is))->Info(f);
 248  }
 249
 250  AliTrack* t=this;
 251  Dump(t,1,f); // Information of all decay products
 252 }
 253 //////////////////////////////////////////////////////////////////////////
 254 void AliTrack::Dump(AliTrack* t,Int_t n,TString f)
 255 {
 256 // Recursively provide the info of all decay levels of this track
 257  AliTrack* td;
 258  for (Int_t id=1; id<=t->GetNdecay(); id++)
 259  {
 260   td=t->GetDecayTrack(id);
 261   if (td)
 262   {
 263    cout << "  ---Level " << n << " sec. track no. " << id << endl;
 264    td->Info(f);
 265    for (Int_t is=1; is<=td->GetNsignals(); is++)
 266    {
 267     ((AliSignal*)td->GetSignal(is))->Info(f);
 268    }
 269
 270    // Go for next decay level of this decay track recursively
 271    Dump(td,n+1,f);
 272   }
 273   else
 274   {
 275    cout << " *AliTrack::Dump* Error : No decay track present." << endl;
 276   }
 277  }
 278 }
 279 //////////////////////////////////////////////////////////////////////////
 280 Double_t AliTrack::GetMomentum()
 281 {
 282 // Provide the value of the track 3-momentum.
 283 // The error can be obtained by invoking GetResultError() after
 284 // invokation of GetMomentum().
 285
 286 // Ali3Vector p=Get3Vector();
 287 // return sqrt(p.Dot(p));
 288  Double_t norm=fV.GetNorm();
 289  return norm;
 290 }
 291 ///////////////////////////////////////////////////////////////////////////
 292 Ali3Vector AliTrack::Get3Momentum()
 293 {
 294 // Provide the track 3-momentum
 295  return (Ali3Vector)Get3Vector();
 296 }
 297 ///////////////////////////////////////////////////////////////////////////
 298 Double_t AliTrack::GetMass()
 299 {
 300 // Provide the particle mass.
 301 // The error can be obtained by invoking GetResultError() after
 302 // invokation of GetMass().
 303  Double_t inv=GetInvariant();
 304  Double_t dinv=GetResultError();
 305  Double_t dm=0;
 306  if (inv >= 0)
 307  {
 308  Double_t m=sqrt(inv);
 309  if (m) dm=dinv/(2.*m);
 310  fDresult=dm;
 311  return m;
 312  }
 313  else
 314  {
 315   cout << "*AliTrack::GetMass* Unphysical situation m**2 = " << inv << endl;
 316   cout << " Value 0 will be returned." << endl;
 317   fDresult=dm;
 318   return 0;
 319  }
 320 }
 321 ///////////////////////////////////////////////////////////////////////////
 322 Float_t AliTrack::GetCharge()
 323 {
 324 // Provide the particle charge
 325  return fQ;
 326 }
 327 ///////////////////////////////////////////////////////////////////////////
 328 Double_t AliTrack::GetEnergy()
 329 {
 330 // Provide the particle's energy.
 331 // The error can be obtained by invoking GetResultError() after
 332 // invokation of GetEnergy().
 333  Double_t E=GetScalar();
 334  if (E>0)
 335  {
 336   return E;
 337  }
 338  else
 339  {
 340   cout << "*AliTrack::GetEnergy* Unphysical situation E = " << E << endl;
 341   cout << " Value 0 will be returned." << endl;
 342   return 0;
 343  }
 344 }
 345 ///////////////////////////////////////////////////////////////////////////
 346 void AliTrack::Decay(Double_t m1,Double_t m2,Double_t thcms,Double_t phicms)
 347 {
 348 // Perform 2-body decay of current track
 349 // m1     : mass of decay product 1
 350 // m2     : mass of decay product 2
 351 // thcms  : cms theta decay angle (in rad.) of m1
 352 // phicms : cms phi decay angle (in rad.) of m1
 353
 354  fNdec=2; // it's a 2-body decay
 355
 356  Double_t M=GetMass();
 357
 358 // Compute the 4-momenta of the decay products in the cms
 359 // Note : p2=p1=pnorm for a 2-body decay
 360  Double_t e1=0;
 361  if (M) e1=((M*M)+(m1*m1)-(m2*m2))/(2.*M);
 362  Double_t e2=0;
 363  if (M) e2=((M*M)+(m2*m2)-(m1*m1))/(2.*M);
 364  Double_t pnorm=(e1*e1)-(m1*m1);
 365  if (pnorm>0.)
 366  {
 367   pnorm=sqrt(pnorm);
 368  }
 369  else
 370  {
 371   pnorm=0;
 372  }
 373
 374  Double_t a[3];
 375  a[0]=pnorm;
 376  a[1]=thcms;
 377  a[2]=phicms;
 378  Ali3Vector p;
 379  p.SetVector(a,"sph");
 380
 381  Ali4Vector pprim1;
 382  pprim1.SetVector(e1,p);
 383  pprim1.SetInvariant(m1*m1);
 384
 385  Ali4Vector pprim2;
 386  p*=-1;
 387  pprim2.SetVector(e2,p);
 388  pprim2.SetInvariant(m2*m2);
 389
 390  // Determine boost parameters from the parent particle
 391  Double_t E=GetEnergy();
 392  p=Get3Vector();
 393  Ali4Vector pmu;
 394  pmu.SetVector(E,p);
 395
 396  AliBoost q;
 397  q.Set4Momentum(pmu);
 398
 399  Ali4Vector p1=q.Inverse(pprim1); // Boost decay product 1
 400  Ali4Vector p2=q.Inverse(pprim2); // Boost decay product 2
 401
 402  // Enter the boosted data into the decay tracks array
 403  if (fDecays)
 404  {
 405   fDecays->Delete();
 406   delete fDecays;
 407  }
 408  fDecays=new TObjArray();
 409
 410  fDecays->Add(new AliTrack);
 411  ((AliTrack*)fDecays->At(0))->Set4Momentum(p1);
 412  ((AliTrack*)fDecays->At(0))->SetMass(m1);
 413  fDecays->Add(new AliTrack);
 414  ((AliTrack*)fDecays->At(1))->Set4Momentum(p2);
 415  ((AliTrack*)fDecays->At(1))->SetMass(m2);
 416 }
 417 ///////////////////////////////////////////////////////////////////////////
 418 Int_t AliTrack::GetNdecay()
 419 {
 420 // Provide the number of decay produced tracks
 421  return fNdec;
 422 }
 423 ///////////////////////////////////////////////////////////////////////////
 424 AliTrack* AliTrack::GetDecayTrack(Int_t j)
 425 {
 426 // Provide decay produced track number j
 427 // Note : j=1 denotes the first decay track
 428  if ((j >= 1) && (j <= fNdec))
 429  {
 430   return (AliTrack*)fDecays->At(j-1);
 431  }
 432  else
 433  {
 434   cout << " *AliTrack* decay track number : " << j << " out of range." << endl;
 435   cout << " -- Decay track number 1 (if any) returned." << endl;
 436   return (AliTrack*)fDecays->At(0);
 437  }
 438 }
 439 ///////////////////////////////////////////////////////////////////////////
 440 void AliTrack::AddSignal(AliSignal& s)
 441 {
 442 // Relate an AliSignal object to this track.
 443  if (!fSignals) fSignals=new TObjArray();
 444  fNsig++;
 445  fSignals->Add(&s);
 446 }
 447 ///////////////////////////////////////////////////////////////////////////
 448 void AliTrack::RemoveSignal(AliSignal& s)
 449 {
 450 // Remove related AliSignal object to this track.
 451  if (fSignals)
 452  {
 453   AliSignal* test=(AliSignal*)fSignals->Remove(&s);
 454   if (test)
 455   {
 456    fNsig--;
 457    fSignals->Compress();
 458   }
 459  }
 460 }
 461 ///////////////////////////////////////////////////////////////////////////
 462 Int_t AliTrack::GetNsignals()
 463 {
 464 // Provide the number of related AliSignals.
 465  return fNsig;
 466 }
 467 ///////////////////////////////////////////////////////////////////////////
 468 AliSignal* AliTrack::GetSignal(Int_t j)
 469 {
 470 // Provide the related AliSignal number j.
 471 // Note : j=1 denotes the first signal.
 472  if ((j >= 1) && (j <= fNsig))
 473  {
 474   return (AliSignal*)fSignals->At(j-1);
 475  }
 476  else
 477  {
 478   cout << " *AliTrack* signal number : " << j << " out of range." << endl;
 479   cout << " -- Signal number 1 (if any) returned." << endl;
 480   return (AliSignal*)fDecays->At(0);
 481  }
 482 }
 483 ///////////////////////////////////////////////////////////////////////////
 484 void AliTrack::SetBeginPoint(AliPosition p)
 485 {
 486 // Store the position of the track begin-point.
 487  fBegin=p;
 488 }
 489 ///////////////////////////////////////////////////////////////////////////
 490 AliPosition AliTrack::GetBeginPoint()
 491 {
 492 // Provide the position of the track begin-point.
 493  return fBegin;
 494 }
 495 ///////////////////////////////////////////////////////////////////////////
 496 void AliTrack::SetEndPoint(AliPosition p)
 497 {
 498 // Store the position of the track end-point.
 499  fEnd=p;
 500 }
 501 ///////////////////////////////////////////////////////////////////////////
 502 AliPosition AliTrack::GetEndPoint()
 503 {
 504 // Provide the position of the track end-point.
 505  return fEnd;
 506 }
 507 ///////////////////////////////////////////////////////////////////////////
 508 void AliTrack::AddMassHypothesis(Double_t prob,Double_t m,Double_t dm)
 509 {
 510 // Add a mass hypothesis for this current track.
 511 // prob=probalility  m=mass value  dm=error on the mass value.
 512 // The default value for the mass error dm is 0.
 513  if (!fMasses) fMasses=new TArrayD();
 514  if (!fDmasses) fDmasses=new TArrayD();
 515  if (!fPmasses) fPmasses=new TArrayD();
 516
 517  fNmasses++;
 518  fMasses->Set(fNmasses);
 519  fDmasses->Set(fNmasses);
 520  fPmasses->Set(fNmasses);
 521
 522  fMasses->AddAt(m,fNmasses-1);
 523  fDmasses->AddAt(dm,fNmasses-1);
 524  fPmasses->AddAt(prob,fNmasses-1);
 525 }
 526 ///////////////////////////////////////////////////////////////////////////
 527 Int_t AliTrack::GetNMassHypotheses()
 528 {
 529 // Provide the number of mass hypotheses for this track.
 530  return fNmasses;
 531 }
 532 ///////////////////////////////////////////////////////////////////////////
 533 Double_t AliTrack::GetMassHypothesis(Int_t j)
 534 {
 535 // Provide the mass of the jth hypothesis for this track.
 536 // Note : the first hypothesis is indicated by j=1.
 537 // Default : j=0 ==> Hypothesis with highest probability.
 538 // The error on the mass can be obtained by invoking GetResultError()
 539 // after invokation of GetMassHypothesis(j).
 540
 541  Double_t m=0,dm=0,prob=0;
 542
 543  // Check validity of index j
 544  if (j<0 || j>fNmasses)
 545  {
 546   cout << " *AliTrack::GetMassHypothesis* Invalid index j : " << j
 547        << " Number of mass hypotheses : " << fNmasses << endl;
 548   fDresult=0;
 549   return 0;
 550  }
 551
 552  // Select mass hypothesis with highest probability
 553  if (j==0)
 554  {
 555   if (fNmasses)
 556   {
 557    m=fMasses->At(0);
 558    dm=fDmasses->At(0);
 559    prob=fPmasses->At(0);
 560    for (Int_t i=1; i<fNmasses; i++)
 561    {
 562     if (fPmasses->At(i)>prob)
 563     {
 564      m=fMasses->At(i);
 565      dm=fDmasses->At(i);
 566     }
 567    }
 568   }
 569   fDresult=dm;
 570   return m;
 571  }
 572
 573  // Provide data of requested mass hypothesis
 574  m=fMasses->At(j-1);
 575  fDresult=fDmasses->At(j-1);
 576  return m;
 577 }
 578 ///////////////////////////////////////////////////////////////////////////
 579 Double_t AliTrack::GetMassHypothesisProb(Int_t j)
 580 {
 581 // Provide the probability of the jth hypothesis for this track.
 582 // Note : the first hypothesis is indicated by j=1.
 583 // Default : j=0 ==> Hypothesis with highest probability.
 584
 585  Double_t prob=0;
 586
 587  // Check validity of index j
 588  if (j<0 || j>fNmasses)
 589  {
 590   cout << " *AliTrack::GetMassHypothesisProb* Invalid index j : " << j
 591        << " Number of mass hypotheses : " << fNmasses << endl;
 592   return 0;
 593  }
 594
 595  // Select mass hypothesis with highest probability
 596  if (j==0)
 597  {
 598   if (fNmasses)
 599   {
 600    prob=fPmasses->At(0);
 601    for (Int_t i=1; i<fNmasses; i++)
 602    {
 603     if (fPmasses->At(i)>prob) prob=fPmasses->At(i);
 604    }
 605   }
 606   return prob;
 607  }
 608
 609  // Provide probability of requested mass hypothesis
 610  prob=fPmasses->At(j-1);
 611  return prob;
 612 }
 613 ///////////////////////////////////////////////////////////////////////////
 614 void AliTrack::SetMass()
 615 {
 616 // Set the mass and error to the value of the hypothesis with highest prob.
 617
 618  Double_t m=0,dm=0,prob=0;
 619
 620  // Select mass hypothesis with highest probability
 621  if (fNmasses)
 622  {
 623   m=fMasses->At(0);
 624   dm=fDmasses->At(0);
 625   prob=fPmasses->At(0);
 626   for (Int_t i=1; i<fNmasses; i++)
 627   {
 628    if (fPmasses->At(i)>prob)
 629    {
 630     m=fMasses->At(i);
 631     dm=fDmasses->At(i);
 632    }
 633   }
 634   SetMass(m,dm);
 635  }
 636  else
 637  {
 638   cout << " *AliTrack::SetMass()* No hypothesis present => No action." << endl;
 639  }
 640 }
 641 ///////////////////////////////////////////////////////////////////////////
 642 void AliTrack::RemoveMassHypothesis(Int_t j)
 643 {
 644 // Remove the jth mass hypothesis for this track.
 645 // Note : the first hypothesis is indicated by j=1.
 646
 647  if (j<=0 || j>fNmasses) // Check validity of index j
 648  {
 649   cout << " *AliTrack::RemoveMassHypothesis* Invalid index j : " << j
 650        << " Number of mass hypotheses : " << fNmasses << endl;
 651  }
 652  else
 653  {
 654   if (j != fNmasses)
 655   {
 656    fMasses->AddAt(fMasses->At(fNmasses-1),j-1);
 657    fDmasses->AddAt(fDmasses->At(fNmasses-1),j-1);
 658    fPmasses->AddAt(fPmasses->At(fNmasses-1),j-1);
 659   }
 660   fMasses->AddAt(0,fNmasses-1);
 661   fDmasses->AddAt(0,fNmasses-1);
 662   fPmasses->AddAt(0,fNmasses-1);
 663   fNmasses--;
 664   fMasses->Set(fNmasses);
 665   fDmasses->Set(fNmasses);
 666   fPmasses->Set(fNmasses);
 667  }
 668 }
 669 ///////////////////////////////////////////////////////////////////////////