[u/mrichter/AliRoot.git] / STEER / AliAODRecoDecay.cxx

/**************************************************************************
 * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/////////////////////////////////////////////////////////////
//
// Base class for AOD reconstructed decay
//
// Author: A.Dainese, andrea.dainese@lnl.infn.it
/////////////////////////////////////////////////////////////

#include <TDatabasePDG.h>
#include <TVector3.h>
#include "AliVParticle.h"
#include "AliAODRecoDecay.h"

ClassImp(AliAODRecoDecay)

//--------------------------------------------------------------------------
AliAODRecoDecay::AliAODRecoDecay() :
  AliVParticle(),
  fSecondaryVtx(0x0),
  fOwnSecondaryVtx(0x0),
  fCharge(0),
  fNProngs(0), fNDCA(0), fNPID(0),
  fPx(0x0), fPy(0x0), fPz(0x0),
  fd0(0x0),
  fDCA(0x0),
  fPID(0x0), 
  fEventNumber(-1),fRunNumber(-1)
{
  //
  // Default Constructor
  //
}
//--------------------------------------------------------------------------
AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
				 Short_t charge,
				 Double_t *px,Double_t *py,Double_t *pz,
				 Double_t *d0) :
  AliVParticle(),
  fSecondaryVtx(vtx2),
  fOwnSecondaryVtx(0x0),
  fCharge(charge),
  fNProngs(nprongs), fNDCA(0), fNPID(0),
  fPx(0x0), fPy(0x0), fPz(0x0),
  fd0(0x0),
  fDCA(0x0),
  fPID(0x0), 
  fEventNumber(-1),fRunNumber(-1)
{
  //
  // Constructor with AliAODVertex for decay vertex
  //

  fPx = new Double_t[GetNProngs()];
  fPy = new Double_t[GetNProngs()];
  fPz = new Double_t[GetNProngs()];
  fd0 = new Double_t[GetNProngs()];
  for(Int_t i=0; i<GetNProngs(); i++) {
    fPx[i] = px[i];
    fPy[i] = py[i];
    fPz[i] = pz[i];
    fd0[i] = d0[i];
  }
}
//--------------------------------------------------------------------------
AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
				 Short_t charge,
				 Double_t *d0) :
  AliVParticle(),
  fSecondaryVtx(vtx2),
  fOwnSecondaryVtx(0x0),
  fCharge(charge),
  fNProngs(nprongs), fNDCA(0), fNPID(0),
  fPx(0x0), fPy(0x0), fPz(0x0),
  fd0(0x0),
  fDCA(0x0),
  fPID(0x0), 
  fEventNumber(-1),fRunNumber(-1)
{
  //
  // Constructor with AliAODVertex for decay vertex and without prongs momenta
  //

  fd0 = new Double_t[GetNProngs()];
  for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
}
//--------------------------------------------------------------------------
AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
  AliVParticle(source),
  fSecondaryVtx(source.fSecondaryVtx),
  fOwnSecondaryVtx(source.fOwnSecondaryVtx),
  fCharge(source.fCharge),
  fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
  fPx(0x0), fPy(0x0), fPz(0x0),
  fd0(0x0), 
  fDCA(0x0),
  fPID(0x0), 
  fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
{
  //
  // Copy constructor
  //
  if(source.GetNProngs()>0) {
    fd0 = new Double32_t[GetNProngs()];
    memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
    if(source.fPx) {
      fPx = new Double32_t[GetNProngs()];
      fPy = new Double32_t[GetNProngs()];
      fPz = new Double32_t[GetNProngs()];
      memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
      memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
      memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
    }
    if(source.fPID) {
      fPID = new Double32_t[fNPID];
      memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
    }
    if(source.fDCA) {
      fDCA = new Double32_t[fNDCA];
      memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
    }
  }
}
//--------------------------------------------------------------------------
AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
{
  //
  // assignment operator
  //
  if(&source == this) return *this;
  fSecondaryVtx = source.fSecondaryVtx;
  fOwnSecondaryVtx = source.fOwnSecondaryVtx;
  fCharge = source.fCharge;
  fNProngs = source.fNProngs;
  fNDCA = source.fNDCA;
  fNPID = source.fNPID;
  fEventNumber = source.fEventNumber;
  fRunNumber = source.fRunNumber;
  if(source.GetNProngs()>0) {
    if(fd0)delete [] fd0; 
    fd0 = new Double32_t[GetNProngs()];
    memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
    if(source.fPx) {
      if(fPx) delete [] fPx; 
      fPx = new Double32_t[GetNProngs()];
      if(fPy) delete [] fPy; 
      fPy = new Double32_t[GetNProngs()];
      if(fPz) delete [] fPz; 
      fPz = new Double32_t[GetNProngs()];
      memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
      memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
      memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
    }
    if(source.fPID) {
      if(fPID) delete [] fPID; 
      fPID = new Double32_t[fNPID];
      memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
    }
    if(source.fDCA) {
      if(fDCA) delete [] fDCA; 
      fDCA = new Double32_t[fNDCA];
      memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
    }
  }
  return *this;
}
//--------------------------------------------------------------------------
AliAODRecoDecay::~AliAODRecoDecay() {
  //  
  // Default Destructor
  //
  if(fPx) { delete [] fPx; fPx=NULL; } 
  if(fPy) { delete [] fPy; fPy=NULL; }
  if(fPz) { delete [] fPz; fPz=NULL; }
  if(fd0) { delete [] fd0; fd0=NULL; }
  if(fPID) { delete [] fPID; fPID=NULL; }
  if(fDCA) { delete [] fDCA; fDCA=NULL; }
  if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::Alpha() const 
{
  //
  // Armenteros-Podolanski alpha for 2-prong decays
  //
  if(GetNProngs()!=2) {
    printf("Can be called only for 2-prong decays");
    return (Double_t)-99999.;
  }
  return 1.-2./(1.+QlProng(0)/QlProng(1));
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const 
{
  //
  // Decay length assuming it is produced at "point" [cm]
  //
  return TMath::Sqrt((point[0]-GetSecVtxX())
		    *(point[0]-GetSecVtxX())
		    +(point[1]-GetSecVtxY())
		    *(point[1]-GetSecVtxY())
		    +(point[2]-GetSecVtxZ())
		    *(point[2]-GetSecVtxZ()));  
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const 
{
  //
  // Decay length in XY assuming it is produced at "point" [cm]
  //
  return TMath::Sqrt((point[0]-GetSecVtxX())
		    *(point[0]-GetSecVtxX())
		    +(point[1]-GetSecVtxY())
		    *(point[1]-GetSecVtxY()));  
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const 
{
  //
  // Cosine of pointing angle in space assuming it is produced at "point"
  //
  TVector3 mom(Px(),Py(),Pz());
  TVector3 fline(GetSecVtxX()-point[0],
		 GetSecVtxY()-point[1],
		 GetSecVtxZ()-point[2]);

  Double_t pta = mom.Angle(fline);

  return TMath::Cos(pta); 
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const 
{
  //
  // Cosine of pointing angle in transverse plane assuming it is produced 
  // at "point"
  //
  TVector3 momXY(Px(),Py(),0.);
  TVector3 flineXY(GetSecVtxX()-point[0],
		   GetSecVtxY()-point[1],
		   0.);

  Double_t ptaXY = momXY.Angle(flineXY);

  return TMath::Cos(ptaXY); 
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const 
{
  //
  // Only for 2-prong decays: 
  // Cosine of decay angle (theta*) in the rest frame of the mother particle
  // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
  // pdgprong0 for prong 0 and pdgprong1 for prong1
  //
  if(GetNProngs()!=2) {
    printf("Can be called only for 2-prong decays");
    return (Double_t)-99999.;
  }
  Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
  Double_t massp[2];
  massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
  massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();

  Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);

  Double_t beta = P()/E(pdgvtx);
  Double_t gamma = E(pdgvtx)/massvtx;

  Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;

  return cts;
}
//---------------------------------------------------------------------------
Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
{
  //
  // Decay time * c assuming it is produced at "point" [cm]
  //
  Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
  return DecayLength(point)*mass/P();
}
//---------------------------------------------------------------------------
Double_t AliAODRecoDecay::E(UInt_t pdg) const 
{
  //
  // Energy
  //
  Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
  return TMath::Sqrt(mass*mass+P()*P());
}
//---------------------------------------------------------------------------
Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const 
{
  //
  // Energy of ip-th prong 
  //
  Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
  return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
}
//---------------------------------------------------------------------------
/*Int_t AliAODRecoDecay::GetIndexProng(Int_t ip) const
{ 
  //
  // Index of prong ip
  //
  if(!GetNProngs()) return 999999;
UShort_t *indices = GetSecondaryVtx()->GetIndices();
  return indices[ip];
}*/
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const 
{
  //
  // Impact parameter in the bending plane of the particle 
  // w.r.t. to "point"
  //
  Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
  k /= Pt()*Pt();
  Double_t dx = GetSecVtxX()-point[0]+k*Px();
  Double_t dy = GetSecVtxY()-point[1]+k*Py();
  Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
  TVector3 mom(Px(),Py(),Pz());
  TVector3 fline(GetSecVtxX()-point[0],
		 GetSecVtxY()-point[1],
		 GetSecVtxZ()-point[2]);
  TVector3 cross = mom.Cross(fline);
  return (cross.Z()>0. ? absImpPar : -absImpPar);
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const 
{
  //
  // Invariant mass for prongs mass hypotheses in pdg array
  //
  if(GetNProngs()!=npdg) {
    printf("npdg != GetNProngs()");
    return (Double_t)-99999.;
  }
  Double_t energysum = 0.;

  for(Int_t i=0; i<GetNProngs(); i++) {
    energysum += EProng(i,pdg[i]);
  }

  Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());

  return mass;
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
				      UInt_t pdg1,UInt_t pdg2) const 
{
  //
  // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
  //
  Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
  Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
                  +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
                  +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
  Double_t mass = TMath::Sqrt(energysum*energysum-psum2);

  return mass;
}
//---------------------------------------------------------------------------
void AliAODRecoDecay::Print(Option_t* /*option*/) const 
{
  //
  // Print some information
  //
  printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
  printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());

  return;
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const 
{
  //
  // Relative angle between two prongs
  //
  TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
  TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));

  Double_t angle = momA.Angle(momB);

  return angle; 
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::QlProng(Int_t ip) const 
{
  //
  // Longitudinal momentum of prong w.r.t. to total momentum
  //
  TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
  TVector3 momTot(Px(),Py(),Pz());

  return mom.Dot(momTot)/momTot.Mag();
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::QtProng(Int_t ip) const 
{
  //
  // Transverse momentum of prong w.r.t. to total momentum  
  //
  TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
  TVector3 momTot(Px(),Py(),Pz());

  return mom.Perp(momTot);
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const 
{
  //
  // Longitudinal momentum of prong w.r.t. to flight line between "point"
  // and fSecondaryVtx
  //
  TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
  TVector3 fline(GetSecVtxX()-point[0],
		 GetSecVtxY()-point[1],
		 GetSecVtxZ()-point[2]);

  return mom.Dot(fline)/fline.Mag();
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const 
{
  //
  // Transverse momentum of prong w.r.t. to flight line between "point" and 
  // fSecondaryVtx 
  //
  TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
  TVector3 fline(GetSecVtxX()-point[0],
		 GetSecVtxY()-point[1],
		 GetSecVtxZ()-point[2]);

  return mom.Perp(fline);
}
//--------------------------------------------------------------------------
Commit	Line	Data
7de7497b	1	/**************************************************************************
	2	* Copyright(c) 1998-2006, 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	// Base class for AOD reconstructed decay
	19	//
	20	// Author: A.Dainese, andrea.dainese@lnl.infn.it
	21	/////////////////////////////////////////////////////////////
	22
	23	#include <TDatabasePDG.h>
	24	#include <TVector3.h>
ff7c57dd	25	#include "AliVParticle.h"
7de7497b	26	#include "AliAODRecoDecay.h"
	27
	28	ClassImp(AliAODRecoDecay)
	29
	30	//--------------------------------------------------------------------------
	31	AliAODRecoDecay::AliAODRecoDecay() :
ff7c57dd	32	AliVParticle(),
7de7497b	33	fSecondaryVtx(0x0),
a11de4a0	34	fOwnSecondaryVtx(0x0),
7de7497b	35	fCharge(0),
	36	fNProngs(0), fNDCA(0), fNPID(0),
	37	fPx(0x0), fPy(0x0), fPz(0x0),
	38	fd0(0x0),
	39	fDCA(0x0),
	40	fPID(0x0),
	41	fEventNumber(-1),fRunNumber(-1)
	42	{
	43	//
	44	// Default Constructor
	45	//
7de7497b	46	}
	47	//--------------------------------------------------------------------------
	48	AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
	49	Short_t charge,
	50	Double_t px,Double_t py,Double_t *pz,
	51	Double_t *d0) :
ff7c57dd	52	AliVParticle(),
7de7497b	53	fSecondaryVtx(vtx2),
a11de4a0	54	fOwnSecondaryVtx(0x0),
7de7497b	55	fCharge(charge),
	56	fNProngs(nprongs), fNDCA(0), fNPID(0),
	57	fPx(0x0), fPy(0x0), fPz(0x0),
	58	fd0(0x0),
	59	fDCA(0x0),
	60	fPID(0x0),
	61	fEventNumber(-1),fRunNumber(-1)
	62	{
	63	//
	64	// Constructor with AliAODVertex for decay vertex
	65	//
	66
	67	fPx = new Double_t[GetNProngs()];
	68	fPy = new Double_t[GetNProngs()];
	69	fPz = new Double_t[GetNProngs()];
	70	fd0 = new Double_t[GetNProngs()];
	71	for(Int_t i=0; i<GetNProngs(); i++) {
	72	fPx[i] = px[i];
	73	fPy[i] = py[i];
	74	fPz[i] = pz[i];
	75	fd0[i] = d0[i];
	76	}
	77	}
	78	//--------------------------------------------------------------------------
	79	AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
	80	Short_t charge,
	81	Double_t *d0) :
ff7c57dd	82	AliVParticle(),
7de7497b	83	fSecondaryVtx(vtx2),
a11de4a0	84	fOwnSecondaryVtx(0x0),
7de7497b	85	fCharge(charge),
	86	fNProngs(nprongs), fNDCA(0), fNPID(0),
	87	fPx(0x0), fPy(0x0), fPz(0x0),
	88	fd0(0x0),
	89	fDCA(0x0),
	90	fPID(0x0),
	91	fEventNumber(-1),fRunNumber(-1)
	92	{
	93	//
	94	// Constructor with AliAODVertex for decay vertex and without prongs momenta
	95	//
	96
	97	fd0 = new Double_t[GetNProngs()];
	98	for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
	99	}
	100	//--------------------------------------------------------------------------
	101	AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
ff7c57dd	102	AliVParticle(source),
7de7497b	103	fSecondaryVtx(source.fSecondaryVtx),
a11de4a0	104	fOwnSecondaryVtx(source.fOwnSecondaryVtx),
7de7497b	105	fCharge(source.fCharge),
	106	fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
	107	fPx(0x0), fPy(0x0), fPz(0x0),
	108	fd0(0x0),
	109	fDCA(0x0),
	110	fPID(0x0),
	111	fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
	112	{
	113	//
	114	// Copy constructor
	115	//
	116	if(source.GetNProngs()>0) {
58b0186f	117	fd0 = new Double32_t[GetNProngs()];
58b0186f	118	memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
7de7497b	119	if(source.fPx) {
58b0186f	120	fPx = new Double32_t[GetNProngs()];
	121	fPy = new Double32_t[GetNProngs()];
	122	fPz = new Double32_t[GetNProngs()];
	123	memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
	124	memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
	125	memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
7de7497b	126	}
7de7497b	127	if(source.fPID) {
58b0186f	128	fPID = new Double32_t[fNPID];
58b0186f	129	memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
7de7497b	130	}
7de7497b	131	if(source.fDCA) {
58b0186f	132	fDCA = new Double32_t[fNDCA];
58b0186f	133	memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
7de7497b	134	}
	135	}
	136	}
	137	//--------------------------------------------------------------------------
	138	AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
	139	{
	140	//
	141	// assignment operator
	142	//
	143	if(&source == this) return *this;
	144	fSecondaryVtx = source.fSecondaryVtx;
a11de4a0	145	fOwnSecondaryVtx = source.fOwnSecondaryVtx;
7de7497b	146	fCharge = source.fCharge;
	147	fNProngs = source.fNProngs;
	148	fNDCA = source.fNDCA;
	149	fNPID = source.fNPID;
	150	fEventNumber = source.fEventNumber;
	151	fRunNumber = source.fRunNumber;
	152	if(source.GetNProngs()>0) {
341952b8	153	if(fd0)delete [] fd0;
58b0186f	154	fd0 = new Double32_t[GetNProngs()];
58b0186f	155	memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
7de7497b	156	if(source.fPx) {
341952b8	157	if(fPx) delete [] fPx;
58b0186f	158	fPx = new Double32_t[GetNProngs()];
341952b8	159	if(fPy) delete [] fPy;
58b0186f	160	fPy = new Double32_t[GetNProngs()];
341952b8	161	if(fPz) delete [] fPz;
58b0186f	162	fPz = new Double32_t[GetNProngs()];
	163	memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
	164	memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
	165	memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
7de7497b	166	}
7de7497b	167	if(source.fPID) {
341952b8	168	if(fPID) delete [] fPID;
58b0186f	169	fPID = new Double32_t[fNPID];
58b0186f	170	memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
7de7497b	171	}
7de7497b	172	if(source.fDCA) {
341952b8	173	if(fDCA) delete [] fDCA;
58b0186f	174	fDCA = new Double32_t[fNDCA];
58b0186f	175	memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
7de7497b	176	}
	177	}
	178	return *this;
	179	}
	180	//--------------------------------------------------------------------------
	181	AliAODRecoDecay::~AliAODRecoDecay() {
	182	//
	183	// Default Destructor
	184	//
a11de4a0	185	if(fPx) { delete [] fPx; fPx=NULL; }
	186	if(fPy) { delete [] fPy; fPy=NULL; }
	187	if(fPz) { delete [] fPz; fPz=NULL; }
	188	if(fd0) { delete [] fd0; fd0=NULL; }
	189	if(fPID) { delete [] fPID; fPID=NULL; }
	190	if(fDCA) { delete [] fDCA; fDCA=NULL; }
	191	if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
7de7497b	192	}
	193	//----------------------------------------------------------------------------
	194	Double_t AliAODRecoDecay::Alpha() const
	195	{
	196	//
	197	// Armenteros-Podolanski alpha for 2-prong decays
	198	//
	199	if(GetNProngs()!=2) {
	200	printf("Can be called only for 2-prong decays");
	201	return (Double_t)-99999.;
	202	}
	203	return 1.-2./(1.+QlProng(0)/QlProng(1));
	204	}
	205	//----------------------------------------------------------------------------
	206	Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
	207	{
	208	//
	209	// Decay length assuming it is produced at "point" [cm]
	210	//
	211	return TMath::Sqrt((point[0]-GetSecVtxX())
	212	*(point[0]-GetSecVtxX())
	213	+(point[1]-GetSecVtxY())
	214	*(point[1]-GetSecVtxY())
	215	+(point[2]-GetSecVtxZ())
	216	*(point[2]-GetSecVtxZ()));
	217	}
	218	//----------------------------------------------------------------------------
	219	Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
	220	{
	221	//
	222	// Decay length in XY assuming it is produced at "point" [cm]
	223	//
	224	return TMath::Sqrt((point[0]-GetSecVtxX())
	225	*(point[0]-GetSecVtxX())
	226	+(point[1]-GetSecVtxY())
	227	*(point[1]-GetSecVtxY()));
	228	}
	229	//----------------------------------------------------------------------------
	230	Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
	231	{
	232	//
	233	// Cosine of pointing angle in space assuming it is produced at "point"
	234	//
	235	TVector3 mom(Px(),Py(),Pz());
	236	TVector3 fline(GetSecVtxX()-point[0],
	237	GetSecVtxY()-point[1],
	238	GetSecVtxZ()-point[2]);
	239
	240	Double_t pta = mom.Angle(fline);
	241
	242	return TMath::Cos(pta);
	243	}
	244	//----------------------------------------------------------------------------
	245	Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
	246	{
	247	//
	248	// Cosine of pointing angle in transverse plane assuming it is produced
	249	// at "point"
	250	//
	251	TVector3 momXY(Px(),Py(),0.);
	252	TVector3 flineXY(GetSecVtxX()-point[0],
	253	GetSecVtxY()-point[1],
	254	0.);
	255
256	Double_t ptaXY = momXY.Angle(flineXY);
257
258	return TMath::Cos(ptaXY);
259	}
260	//----------------------------------------------------------------------------
261	Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
262	{
263	//
264	// Only for 2-prong decays:
265	// Cosine of decay angle (theta*) in the rest frame of the mother particle
266	// for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
267	// pdgprong0 for prong 0 and pdgprong1 for prong1
268	//
269	if(GetNProngs()!=2) {
270	printf("Can be called only for 2-prong decays");
271	return (Double_t)-99999.;
272	}
273	Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
274	Double_t massp[2];
275	massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
276	massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
277
278	Double_t pStar = TMath::Sqrt(TMath::Power(massvtxmassvtx-massp[0]massp[0]-massp[1]massp[1],2.)-4.massp[0]massp[0]massp[1]massp[1])/(2.massvtx);
279
280	Double_t beta = P()/E(pdgvtx);
281	Double_t gamma = E(pdgvtx)/massvtx;
282
283	Double_t cts = (QlProng(ip)/gamma-betaTMath::Sqrt(pStarpStar+massp[ip]*massp[ip]))/pStar;
284
285	return cts;
286	}
287	//---------------------------------------------------------------------------
288	Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
289	{
290	//
291	// Decay time * c assuming it is produced at "point" [cm]
292	//
293	Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
294	return DecayLength(point)*mass/P();
295	}
296	//---------------------------------------------------------------------------
297	Double_t AliAODRecoDecay::E(UInt_t pdg) const
298	{
299	//
300	// Energy
301	//
302	Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
303	return TMath::Sqrt(massmass+P()P());
304	}
305	//---------------------------------------------------------------------------
306	Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
307	{
308	//
309	// Energy of ip-th prong
310	//
311	Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
312	return TMath::Sqrt(massmass+PProng(ip)PProng(ip));
313	}
314	//---------------------------------------------------------------------------
315	/*Int_t AliAODRecoDecay::GetIndexProng(Int_t ip) const
316	{
317	//
318	// Index of prong ip
319	//
320	if(!GetNProngs()) return 999999;
5cc73331	321	UShort_t *indices = GetSecondaryVtx()->GetIndices();
7de7497b	322	return indices[ip];
	323	}*/
	324	//----------------------------------------------------------------------------
	325	Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
	326	{
	327	//
	328	// Impact parameter in the bending plane of the particle
	329	// w.r.t. to "point"
	330	//
	331	Double_t k = -(GetSecVtxX()-point[0])Px()-(GetSecVtxY()-point[1])Py();
	332	k /= Pt()*Pt();
	333	Double_t dx = GetSecVtxX()-point[0]+k*Px();
	334	Double_t dy = GetSecVtxY()-point[1]+k*Py();
	335	Double_t absImpPar = TMath::Sqrt(dxdx+dydy);
	336	TVector3 mom(Px(),Py(),Pz());
	337	TVector3 fline(GetSecVtxX()-point[0],
	338	GetSecVtxY()-point[1],
	339	GetSecVtxZ()-point[2]);
	340	TVector3 cross = mom.Cross(fline);
	341	return (cross.Z()>0. ? absImpPar : -absImpPar);
	342	}
	343	//----------------------------------------------------------------------------
	344	Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
	345	{
	346	//
	347	// Invariant mass for prongs mass hypotheses in pdg array
	348	//
	349	if(GetNProngs()!=npdg) {
	350	printf("npdg != GetNProngs()");
	351	return (Double_t)-99999.;
	352	}
	353	Double_t energysum = 0.;
	354
	355	for(Int_t i=0; i<GetNProngs(); i++) {
	356	energysum += EProng(i,pdg[i]);
	357	}
	358
	359	Double_t mass = TMath::Sqrt(energysumenergysum-P()P());
	360
	361	return mass;
	362	}
	363	//----------------------------------------------------------------------------
	364	Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
	365	UInt_t pdg1,UInt_t pdg2) const
	366	{
	367	//
	368	// 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
	369	//
	370	Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
	371	Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
	372	+(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
	373	+(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
	374	Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
	375
	376	return mass;
	377	}
	378	//---------------------------------------------------------------------------
	379	void AliAODRecoDecay::Print(Option_t* /option/) const
	380	{
	381	//
	382	// Print some information
	383	//
	384	printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
	385	printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
386
387	return;
388	}
389	//----------------------------------------------------------------------------
390	Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
391	{
392	//
393	// Relative angle between two prongs
394	//
395	TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
396	TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
397
398	Double_t angle = momA.Angle(momB);
399
400	return angle;
401	}
402	//----------------------------------------------------------------------------
403	Double_t AliAODRecoDecay::QlProng(Int_t ip) const
404	{
405	//
406	// Longitudinal momentum of prong w.r.t. to total momentum
407	//
408	TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
409	TVector3 momTot(Px(),Py(),Pz());
410
411	return mom.Dot(momTot)/momTot.Mag();
412	}
413	//----------------------------------------------------------------------------
414	Double_t AliAODRecoDecay::QtProng(Int_t ip) const
415	{
416	//
417	// Transverse momentum of prong w.r.t. to total momentum
418	//
419	TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
420	TVector3 momTot(Px(),Py(),Pz());
421
422	return mom.Perp(momTot);
423	}
424	//----------------------------------------------------------------------------
425	Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
426	{
427	//
428	// Longitudinal momentum of prong w.r.t. to flight line between "point"
429	// and fSecondaryVtx
430	//
431	TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
432	TVector3 fline(GetSecVtxX()-point[0],
433	GetSecVtxY()-point[1],
434	GetSecVtxZ()-point[2]);
435
436	return mom.Dot(fline)/fline.Mag();
437	}
438	//----------------------------------------------------------------------------
439	Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
440	{
441	//
442	// Transverse momentum of prong w.r.t. to flight line between "point" and
443	// fSecondaryVtx
444	//
445	TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
446	TVector3 fline(GetSecVtxX()-point[0],
447	GetSecVtxY()-point[1],
448	GetSecVtxZ()-point[2]);
449
450	return mom.Perp(fline);
451	}
452	//--------------------------------------------------------------------------