-/**************************************************************************
- * Copyright(c) 1998-1999, 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. *
- **************************************************************************/
-
-//-------------------------------------------------------------------------
-// Class AliRsnDaughter
-// ----------------------
-// A simple object which describes a reconstructed track
-// with some references to its related generated particle
-// and some facilities which could help in composing its
-// 4-momentum, for resonance study.
-//
-// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
-//-------------------------------------------------------------------------
+//
+// Class AliRsnDaughter
+//
+// Interface to candidate daughters of a resonance (tracks).
+// Points to the source of information, which is generally an AliVParticle-derived object
+// and contains few internal data-members to store "on fly" some important information
+// for the computations required during resonance analysis.
+// ---
+// Since the package revision, this object is not supposed to be stacked in memory
+// but created "on fly" during analysis and used just for computations, as an interface.
+//
+// authors: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
+// M. Vala (martin.vala@cern.ch)
+//
+#include <Riostream.h>
#include <TParticle.h>
-#include <TString.h>
-#include "AliLog.h"
+#include "AliStack.h"
#include "AliESDtrack.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+
+#include "AliRsnPIDDefESD.h"
#include "AliRsnDaughter.h"
ClassImp(AliRsnDaughter)
-//--------------------------------------------------------------------------------------------------------
-AliRsnDaughter::AliRsnDaughter() :
- TObject(),
- fSign((Char_t)0),
- fPDG((UShort_t)0),
- fIndex((UShort_t)0),
- fMass(0.0),
- fLabel(-1),
- fTruePDG((Short_t)0),
- fMother(-1),
- fMotherPDG((Short_t)0)
-{
-//
+AliRsnDaughter::EPIDMethod AliRsnDaughter::fgPIDMethod = AliRsnDaughter::kRealistic;
+
+//_____________________________________________________________________________
+AliRsnDaughter::AliRsnDaughter(AliVParticle *ref, TParticle *refMC) :
+ fOK((ref != 0)),
+ fKinkIndex(0),
+ fParticle(refMC),
+ fMotherPDG(0),
+ fStatus(0),
+ fDr(0.0),
+ fDz(0.0),
+ fReqPID(AliPID::kUnknown),
+ fRef(ref)
+{
+//
// Default constructor.
-// Its unique argument defines how many PID weights are allowed.
-// Actually, it should be the same as AliESDtrack::kSPECIES (=5).
-//
-
- Int_t i;
- for (i = 0; i < AliPID::kSPECIES; i++) {
- if (i < 3) {
- fP[i] = 0.0;
- fV[i] = 0.0;
- }
- fPIDwgt[i] = 0.0;
- }
-}
-//--------------------------------------------------------------------------------------------------------
-AliRsnDaughter::AliRsnDaughter(Int_t label, UShort_t index, Double_t *p, Double_t *v, Char_t sign) :
- TObject(),
- fSign(sign),
- fPDG((UShort_t)0),
- fIndex(index),
- fMass(0.0),
- fLabel(label),
- fTruePDG((Short_t)0),
- fMother(-1),
- fMotherPDG((Short_t)0)
-{
-//
-// [PRIVATE]
-// Constructor with arguments.
-// To create an object with all its standard data members filled (the ones coming from ESD)
-//
- Int_t i;
- for (i = 0; i < AliPID::kSPECIES; i++) {
- if (i < 3) {
- fP[i] = p[i];
- fV[i] = v[i];
- }
- fPIDwgt[i] = 0.0;
- }
-}
-//--------------------------------------------------------------------------------------------------------
-AliRsnDaughter * AliRsnDaughter::Adopt(AliESDtrack* esdTrack, Int_t index)
-{
-//
-// Copies reconstructed data from an AliESDtrack:
-//
-// - charge sign
-// - momentum
-// - point of closest approach to primary vertex
-// - ESD pid weights
-// - track label (AliESDtrack::GetLabel())
-//
-// Given that it happened that tracks have exactly zero momentum and DCA point,
-// these quantities are checked and skipped when caught.
-//
- Double_t p[3], v[3], pid[AliPID::kSPECIES];
-
- esdTrack->GetPxPyPz(p);
- esdTrack->GetXYZ(v);
- esdTrack->GetESDpid(pid);
-
- if (p[0] == 0.0 || p[1] == 0.0 || p[2] == 0.0) return 0x0;
- if (v[0] == 0.0 || v[1] == 0.0 || v[2] == 0.0) return 0x0;
-
- AliRsnDaughter *out = new AliRsnDaughter(esdTrack->GetLabel(), (UShort_t)index, p, v, (Char_t)esdTrack->GetSign());
-
- // store PID weights
- out->SetPIDweights(pid);
-
- return out;
-}
-//--------------------------------------------------------------------------------------------------------
-AliRsnDaughter * AliRsnDaughter::Adopt(TParticle* particle, Int_t label)
-{
-//
-// Copies data from a generated particle:
-//
-// - PDG code
-// - charge sign
-// - momentum
-// - production vertex
-// - GEANT label of mother track
-//
- Int_t pdg;
- Char_t sign;
- Double_t p[3], v[3];
-
- // get particle sign form the sign of PDG code
- pdg = particle->GetPdgCode();
- if (TMath::Abs(pdg) < 20) {
- if (pdg > 0) sign = -1; else sign = 1;
- }
- else if (TMath::Abs(pdg) < 3000) {
- if (pdg > 0) sign = 1; else sign = -1;
- }
- else {
- return 0x0;
- }
-
- p[0] = particle->Px();
- p[1] = particle->Py();
- p[2] = particle->Pz();
-
- v[0] = particle->Vx();
- v[1] = particle->Vy();
- v[2] = particle->Vz();
-
- AliRsnDaughter *out = new AliRsnDaughter(label, (UShort_t)TMath::Abs(label), p, v, sign);
-
- // set simulation data
- out->SetPDG(TMath::Abs(pdg));
- out->SetTruePDG(TMath::Abs(pdg));
- out->SetMother(particle->GetFirstMother());
-
- return out;
-}
-//--------------------------------------------------------------------------------------------------------
-void AliRsnDaughter::Print(Option_t *option) const
-{
-//
-// Print informations about track.
-// All options are used to add some specifical information as required:
-// "S" --> charge sign
-// "L" --> track label
-// "I" --> assigned PID (PDG code)
-// "T" --> true PDG code
-// "V" --> coordinates of track vertex
-// "P" --> coordinates of track momentum
-// "M" --> mother label
-// "N" --> mother PDG code
-// "W" --> ESD PID weights
-//
-
- TString output("Track info: ");
- TString opt(option);
- opt.ToUpper();
-
- if (opt.Contains("S")) {
- output.Append(Form("sign = %d -- ", (Int_t)fSign));
- }
- if (opt.Contains("L")) {
- output.Append(Form("label = %d -- ", fLabel));
- }
- if (opt.Contains("I")) {
- output.Append(Form("PDG = %d -- ", fPDG));
- }
- if (opt.Contains("T")) {
- output.Append(Form("true PDG = %d -- ", fTruePDG));
- }
- if (opt.Contains("V")) {
- output.Append(Form("v = %f, %f, %f -- ", fV[0], fV[1], fV[2]));
- }
- if (opt.Contains("P")) {
- output.Append(Form("p = %f, %f, %f -- ", fP[0], fP[1], fP[2]));
- }
- if (opt.Contains("M")) {
- output.Append(Form("mum = %d -- ", fMother));
- }
- if (opt.Contains("N")) {
- output.Append(Form("mum PDG = %d -- ", fMotherPDG));
- }
- if (opt.Contains("W")) {
- output.Append(Form("PID wgts (e, mu, pi, K, p) = %f, %f, %f, %f, %f", fPIDwgt[0], fPIDwgt[1], fPIDwgt[2], fPIDwgt[3], fPIDwgt[4]));
- }
-
- AliInfo(output.Data());
-}
-//--------------------------------------------------------------------------------------------------------
-AliRsnDaughter AliRsnDaughter::Sum(AliRsnDaughter t1, AliRsnDaughter t2)
-{
-//
-// Builds a new AliRsnDaughter object with the sum of momenta of two particles.
-//
- // create new AliRsnDaughter with default useless values for datamembers
- AliRsnDaughter out;
-
- // if the summed particles are daughters of the same resonance
- // their common mother label becomes the label of the sum
- Int_t mum1 = t1.GetMother();
- Int_t mum2 = t2.GetMother();
- if (mum1 == mum2) {
- out.SetLabel(mum1);
- out.SetMotherPDG(t1.GetMotherPDG());
- }
- else {
- out.SetLabel(-1);
- out.SetMotherPDG(0);
- }
-
- // compute total 4-momentum
- Double_t etot = t1.GetEnergy() + t2.GetEnergy();
- Double_t pxTot = t1.GetPx() + t2.GetPx();
- Double_t pyTot = t1.GetPy() + t2.GetPy();
- Double_t pzTot = t1.GetPz() + t2.GetPz();
- Double_t mass = TMath::Sqrt(etot*etot - pxTot*pxTot - pyTot*pyTot - pzTot*pzTot);
-
- out.SetPxPyPz(pxTot, pyTot, pzTot);
- out.SetMass(mass);
-
- return out;
+//
+}
+
+//_____________________________________________________________________________
+AliRsnDaughter::AliRsnDaughter(const AliRsnDaughter ©) :
+ TObject(copy),
+ fOK(copy.fOK),
+ fKinkIndex(copy.fKinkIndex),
+ fParticle(copy.fParticle),
+ fMotherPDG(copy.fMotherPDG),
+ fStatus(copy.fStatus),
+ fDr(copy.fDr),
+ fDz(copy.fDz),
+ fReqPID(copy.fReqPID),
+ fRef(copy.fRef)
+{
+//
+// Copy constructor.
+// Pointers are NOT duplicated.
+//
+}
+
+//_____________________________________________________________________________
+AliRsnDaughter& AliRsnDaughter::operator=(const AliRsnDaughter ©)
+{
+//
+// Assignment operator.
+//
+
+ (TObject)(*this) = (TObject)copy;
+
+ fOK = copy.fOK;
+ fKinkIndex = copy.fKinkIndex;
+ fParticle = copy.fParticle;
+ fMotherPDG = copy.fMotherPDG;
+ fStatus = copy.fStatus;
+ fDr = copy.fDr;
+ fDz = copy.fDz;
+
+ fReqPID = copy.fReqPID;
+
+ fRef = copy.fRef;
+
+ return (*this);
+}
+
+//_____________________________________________________________________________
+AliRsnDaughter::~AliRsnDaughter()
+{
+//
+// Destructor.
+// Since pointers do not allocate new objects, nothing is done.
+//
+}
+
+//_____________________________________________________________________________
+void AliRsnDaughter::RotateP
+(Double_t angle, Double_t &x, Double_t &y, Bool_t isDegrees) const
+{
+//
+// Rotate the transverse momentum by an angle (in DEGREES)
+// around Z axis (does not change the Z component).
+// Rotated values are stored in the two arguments passed by reference.
+//
+
+ if (isDegrees) angle *= TMath::DegToRad();
+
+ Double_t s = TMath::Sin(angle);
+ Double_t c = TMath::Cos(angle);
+
+ x = c*Px() - s*Py();
+ y = s*Px() + c*Py();
+}
+
+//_____________________________________________________________________________
+Double_t AliRsnDaughter::AngleTo(AliRsnDaughter d, Bool_t outInDegrees)
+{
+//
+// Compute angle between the vector momentum of this
+// and the one of argument.
+//
+
+ Double_t arg, dot, ptot2 = P2() * d.P2();
+
+ if (ptot2 <= 0) {
+ return 0.0;
+ } else {
+ dot = Px()*d.Px() + Py()*d.Py() + Pz()*d.Pz();
+ arg = dot / TMath::Sqrt(ptot2);
+ if (arg > 1.0) arg = 1.0;
+ if (arg < -1.0) arg = -1.0;
+ if (outInDegrees) return TMath::ACos(arg) * TMath::RadToDeg();
+ else return TMath::ACos(arg);
+ }
+}
+
+//_____________________________________________________________________________
+Int_t AliRsnDaughter::GetID() const
+{
+//
+// Return reference index, using the "GetID" method
+// of the possible source object.
+//
+
+ AliESDtrack *esd = dynamic_cast<AliESDtrack*>(fRef);
+ if (esd) return esd->GetID();
+
+ AliAODTrack *aod = dynamic_cast<AliAODTrack*>(fRef);
+ if (aod) return aod->GetID();
+
+ return GetLabel();
+}
+
+//_____________________________________________________________________________
+AliPID::EParticleType AliRsnDaughter::RealisticPID() const
+{
+//
+// Return the "realistic" PID of this track,
+// i.e. the particle species to which corresponds the largest PID probability.
+//
+
+ AliPID::EParticleType pid = AliPID::kElectron;
+ Double_t prob = fPID[0];
+
+ Int_t i;
+ for (i = 1; i < AliPID::kSPECIES; i++) {
+ if (fPID[i] > prob) {
+ prob = fPID[i];
+ pid = (AliPID::EParticleType)i;
+ }
+ }
+
+ return pid;
+}
+
+//_____________________________________________________________________________
+AliPID::EParticleType AliRsnDaughter::PerfectPID() const
+{
+//
+// Return the "perfect" PID of this track,
+// reading it from the MC information, if available.
+//
+
+ if (!fParticle) return AliPID::kUnknown;
+
+ Int_t absPDG = TMath::Abs(fParticle->GetPdgCode());
+ switch (absPDG) {
+ case 11:
+ return AliPID::kElectron;
+ case 13:
+ return AliPID::kMuon;
+ case 211:
+ return AliPID::kPion;
+ case 321:
+ return AliPID::kKaon;
+ case 2212:
+ return AliPID::kProton;
+ default:
+ AliDebug(2, Form("PDG code = %d not recognized. Return 'AliPID::kUnknown'", absPDG));
+ return AliPID::kUnknown;
+ }
+}
+
+//_____________________________________________________________________________
+AliPID::EParticleType AliRsnDaughter::PIDType(Double_t &prob) const
+{
+//
+// Return the PID type according to the selected method
+// in the argument passed by reference, the probability is stored.
+// It will be realistic for realistic PID and 1 for perfect PID.
+//
+
+ AliPID::EParticleType pid = AssignedPID();
+
+ prob = 1.0;
+ if (fgPIDMethod == kRealistic) prob = fPID[(Int_t)pid];
+
+ return pid;
+}
+
+//_____________________________________________________________________________
+AliPID::EParticleType AliRsnDaughter::AssignedPID() const
+{
+//
+// Return the PID type according to the selected method
+// in the argument passed by reference, the probability is stored.
+// It will be realistic for realistic PID and 1 for perfect PID.
+//
+
+ switch (fgPIDMethod) {
+ case kNoPID:
+ return AliPID::kUnknown;
+ case kPerfect:
+ return PerfectPID();
+ case kRealistic:
+ return RealisticPID();
+ default:
+ AliWarning("PID method not properly set. Returning realistic PID");
+ return RealisticPID();
+ }
+}
+
+//_____________________________________________________________________________
+Bool_t AliRsnDaughter::CombineWithPriors(const Double_t *priors, AliRsnPIDDefESD *pidDef)
+{
+//
+// Combine current PID weights (assumed to be them) with prior probs
+//
+
+ Int_t i;
+ Double_t sum = 0.0;
+
+ // get PID weights according to definition
+ // if the reference is not ESD or the pidDef is null
+ // of it is not null but is requires the ESD pid,
+ // the standard PID value is used, otherwise
+ if (pidDef && GetRefESD()) {
+ pidDef->ComputeWeights(GetRefESD(), fPID);
+ } else {
+ if (fRef->PID())
+ for (i = 0; i < AliPID::kSPECIES; i++) fPID[i] = fRef->PID()[i];
+ else if (fParticle) {
+ // if the PID() returns 0x0, this is a MC particle
+ // and we set the weight corresponding to its species to 1
+ // and the others to zero
+ for (i = 0; i < AliPID::kSPECIES; i++) fPID[i] = 0.0;
+ i = (Int_t)AliRsnDaughter::InternalType(TMath::Abs(fParticle->GetPdgCode()));
+ fPID[i] = 1.0;
+ }
+ }
+
+ // multiply weights and priors
+ for (i = 0; i < AliPID::kSPECIES; i++) {
+ fPID[i] *= priors[i];
+ sum += fPID[i];
+ }
+ if (sum <= 0.0) {
+ if (sum < 0.0) AliError(Form("Sum of weights = %f < 0", sum));
+ return kFALSE;
+ }
+
+ // normalize
+ for (i = 0; i < AliPID::kSPECIES; i++) fPID[i] /= sum;
+
+ return kTRUE;
+}
+
+//_____________________________________________________________________________
+void AliRsnDaughter::FindMotherPDG(AliStack *const stack)
+{
+//
+// Searches the stack to find the mother and retrieve its PDG code.
+//
+
+ if (!stack || !fParticle) return;
+
+ const Int_t mLabel = fParticle->GetFirstMother();
+ if (mLabel < 0) {
+ fMotherPDG = 0;
+ } else {
+ TParticle *mum = stack->Particle(mLabel);
+ if (mum) fMotherPDG = mum->GetPdgCode();
+ else fMotherPDG = 0;
+ }
+}
+
+//_____________________________________________________________________________
+Double_t AliRsnDaughter::GetMCEnergy(Double_t mass)
+{
+//
+// Uses the argument to compute 4-momentum energy
+//
+
+ if (!fParticle) return 0.0;
+
+ Double_t p2 = fParticle->Px()*fParticle->Px();
+ p2 += fParticle->Py()*fParticle->Py();
+ p2 += fParticle->Pz()*fParticle->Pz();
+
+ return TMath::Sqrt(mass*mass + p2);
+}
+
+//_____________________________________________________________________________
+void AliRsnDaughter::FindKinkIndex(const AliESDtrack *esdTrack)
+{
+//
+// Assign kink index from an ESD track
+//
+
+ Int_t i, ik[3];
+ for (i = 0; i < 3; i++) ik[i] = esdTrack->GetKinkIndex(i);
+
+ if (ik[0] < 0 || ik[1] < 0 || ik[2] < 0) {
+ SetKinkMother();
+ } else if (ik[0] > 0 || ik[1] > 0 || ik[2] > 0) {
+ SetKinkDaughter();
+ } else SetNoKink();
+}
+
+//_____________________________________________________________________________
+void AliRsnDaughter::FindKinkIndex(AliAODEvent *const event)
+{
+//
+// Assign kink index from an AOD event
+//
+
+ Int_t iv, id, nD, nV = event->GetNumberOfVertices();
+ for (iv = 0; iv < nV; iv++) {
+ AliAODVertex *v = event->GetVertex(iv);
+ AliAODVertex::AODVtx_t type = (AliAODVertex::AODVtx_t)v->GetType();
+ if (type != AliAODVertex::kKink) continue;
+ AliAODTrack *mother = (AliAODTrack*)v->GetParent();
+ if (mother == (AliAODTrack*)fRef) {
+ SetKinkMother();
+ return;
+ } else {
+ nD = v->GetNDaughters();
+ for (id = 0; id < nD; id++) {
+ AliAODTrack *son = (AliAODTrack*)v->GetDaughter(id);
+ if (son == (AliAODTrack*)fRef) {
+ SetKinkDaughter();
+ return;
+ }
+ }
+ }
+ }
+
+ SetNoKink();
+}
+
+//_____________________________________________________________________________
+void AliRsnDaughter::Reset()
+{
+//
+// Reset this track to meaningless values
+//
+
+ fOK = kFALSE;
+ fKinkIndex = 0;
+ fParticle = 0x0;
+ fMotherPDG = 0;
+ fStatus = 0;
+ fRef = 0x0;
+}
+
+//_____________________________________________________________________________
+void AliRsnDaughter::Print(Option_t * const option) const
+{
+//
+// Prints the values of data members, using the options:
+// - P --> momentum
+// - V --> DCA vertex
+// - C --> electric charge
+// - F --> flags
+// - I --> identification (PID, probability and mass)
+// - W --> PID weights
+// - M --> Montecarlo
+// - L --> index & label
+// - A --> angles
+// - ALL --> All oprions switched on
+//
+// Index and label are printed by default.
+//
+
+ TString opt(option);
+ opt.ToUpper();
+
+ if (opt.Contains("L") || opt.Contains("ALL")) {
+ cout << ".......Index : " << GetID() << endl;
+ cout << ".......Label : " << GetLabel() << endl;
+ }
+ if (opt.Contains("P") || opt.Contains("ALL")) {
+ cout << ".......Px, Py, Pz, Pt : " << Px() << ' ' << Py() << ' ' << Pz() << ' ' << Pt() << endl;
+ }
+ if (opt.Contains("A") || opt.Contains("ALL")) {
+ cout << ".......Phi, Theta : " << Phi() << ' ' << Theta() << endl;
+ }
+ if (opt.Contains("V") || opt.Contains("ALL")) {
+ cout << ".......Vx, Vy, Vz : " << Xv() << ' ' << Yv() << ' ' << Zv() << endl;
+ }
+ if (opt.Contains("I") || opt.Contains("ALL")) {
+ AliPID::EParticleType type;
+ Double_t prob;
+ type = PIDType(prob);
+ cout << ".......PID & prob : " << AliPID::ParticleName(type) << ' ' << prob << endl;
+ }
+ if (opt.Contains("C") || opt.Contains("ALL")) {
+ cout << ".......Charge : " << Charge() << endl;
+ }
+ if (opt.Contains("F") || opt.Contains("ALL")) {
+ cout << ".......Flags : " << fStatus << endl;
+ }
+ if (opt.Contains("W") || opt.Contains("ALL")) {
+ cout << ".......Weights : ";
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++) cout << fPID[i] << ' ';
+ cout << endl;
+ }
+ if (opt.Contains("M") || opt.Contains("ALL")) {
+ if (fParticle) {
+ cout << ".......PDG code : " << fParticle->GetPdgCode() << endl;
+ cout << ".......Mother (label) : " << fParticle->GetFirstMother() << endl;
+ cout << ".......Mother (PDG code): " << fMotherPDG << endl;
+ } else {
+ cout << ".......MC info not present" << endl;
+ }
+ }
+}
+
+//_____________________________________________________________________________
+AliPID::EParticleType AliRsnDaughter::InternalType(Int_t pdg)
+{
+//
+// Return the internal enum value corresponding to the PDG
+// code passed as argument, if possible.
+// Otherwise, returns 'AliPID::kSPECIES' by default.
+//
+
+ AliPID::EParticleType value;
+ Int_t absPDG = TMath::Abs(pdg);
+
+ switch (absPDG) {
+ case 11:
+ value = AliPID::kElectron;
+ break;
+ case 13:
+ value = AliPID::kMuon;
+ break;
+ case 211:
+ value = AliPID::kPion;
+ break;
+ case 321:
+ value = AliPID::kKaon;
+ break;
+ case 2212:
+ value = AliPID::kProton;
+ break;
+ default:
+ value = AliPID::kUnknown;
+ }
+ return value;
+}
+
+//_____________________________________________________________________________
+const char* AliRsnDaughter::MethodName(EPIDMethod method)
+{
+//
+// Returns a string with the method name
+//
+
+ switch (method)
+ {
+ case kNoPID:
+ return "No PID";
+ case kRealistic:
+ return "Realistic";
+ case kPerfect:
+ return "Perfect";
+ default:
+ return "Unknown";
+ }
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff