PWG2/RESONANCES/AliRsnDaughter.cxx

   1 //
   2 // Class AliRsnDaughter
   3 //
   4 // Interface to candidate daughters of a resonance (tracks).
   5 // Points to the source of information, which is generally an AliVParticle-derived object
   6 // and contains few internal data-members to store "on fly" some important information
   7 // for the computations required during resonance analysis.
   8 // ---
   9 // Since the package revision, this object is not supposed to be stacked in memory
  10 // but created "on fly" during analysis and used just for computations, as an interface.
  11 //
  12 // authors: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
  13 //          M. Vala (martin.vala@cern.ch)
  14 //
  15
  16 #include <Riostream.h>
  17
  18 #include <TParticle.h>
  19 #include <TString.h>
  20
  21 #include "AliLog.h"
  22 #include "AliStack.h"
  23 #include "AliESDtrack.h"
  24 #include "AliAODEvent.h"
  25 #include "AliAODVertex.h"
  26 #include "AliAODTrack.h"
  27
  28 #include "AliRsnDaughter.h"
  29
  30 ClassImp(AliRsnDaughter)
  31
  32 AliRsnDaughter::EPIDMethod AliRsnDaughter::fgPIDMethod = AliRsnDaughter::kRealistic;
  33
  34 //_____________________________________________________________________________
  35 AliRsnDaughter::AliRsnDaughter(AliVParticle *ref, TParticle *refMC) :
  36   fOK((ref != 0)),
  37   fKinkIndex(0),
  38   fParticle(refMC),
  39   fMotherPDG(0),
  40   fStatus(0),
  41   fRef(ref)
  42 {
  43 //
  44 // Default constructor.
  45 //
  46 }
  47
  48 //_____________________________________________________________________________
  49 AliRsnDaughter::AliRsnDaughter(const AliRsnDaughter &copy) :
  50   TObject(copy),
  51   fOK(copy.fOK),
  52   fKinkIndex(copy.fKinkIndex),
  53   fParticle(copy.fParticle),
  54   fMotherPDG(copy.fMotherPDG),
  55   fStatus(copy.fStatus),
  56   fRef(copy.fRef)
  57 {
  58 //
  59 // Copy constructor.
  60 // Pointers are NOT duplicated.
  61 //
  62 }
  63
  64 //_____________________________________________________________________________
  65 AliRsnDaughter& AliRsnDaughter::operator=(const AliRsnDaughter &copy)
  66 {
  67 //
  68 // Assignment operator.
  69 //
  70
  71   (TObject)(*this) = (TObject)copy;
  72
  73   fOK = copy.fOK;
  74   fKinkIndex = copy.fKinkIndex;
  75   fParticle  = copy.fParticle;
  76   fMotherPDG = copy.fMotherPDG;
  77   fStatus = copy.fStatus;
  78
  79   fRef = copy.fRef;
  80
  81   return (*this);
  82 }
  83
  84 //_____________________________________________________________________________
  85 AliRsnDaughter::~AliRsnDaughter()
  86 {
  87 //
  88 // Destructor.
  89 // Since pointers do not allocate new objects, nothing is done.
  90 //
  91 }
  92
  93 //_____________________________________________________________________________
  94 void AliRsnDaughter::RotateP
  95 (Double_t angle, Double_t &x, Double_t &y, Bool_t isDegrees)
  96 {
  97 //
  98 // Rotate the transverse momentum by an angle (in DEGREES)
  99 // around Z axis (does not change the Z component).
 100 // Rotated values are stored in the two arguments passed by reference.
 101 //
 102
 103   if (isDegrees) angle *= TMath::DegToRad();
 104
 105   Double_t s = TMath::Sin(angle);
 106   Double_t c = TMath::Cos(angle);
 107
 108   x = c*Px() - s*Py();
 109   y = s*Px() + c*Py();
 110 }
 111
 112 //_____________________________________________________________________________
 113 Double_t AliRsnDaughter::AngleTo(AliRsnDaughter d, Bool_t outInDegrees)
 114 {
 115 //
 116 // Compute angle between the vector momentum of this
 117 // and the one of argument.
 118 //
 119
 120   Double_t arg, dot, ptot2 = P2() * d.P2();
 121
 122   if (ptot2 <= 0) {
 123     return 0.0;
 124   } else {
 125     dot = Px()*d.Px() + Py()*d.Py() + Pz()*d.Pz();
 126     arg = dot / TMath::Sqrt(ptot2);
 127     if (arg >  1.0) arg =  1.0;
 128     if (arg < -1.0) arg = -1.0;
 129     if (outInDegrees) return TMath::ACos(arg) * TMath::RadToDeg();
 130     else return TMath::ACos(arg);
 131   }
 132 }
 133
 134 //_____________________________________________________________________________
 135 Int_t AliRsnDaughter::GetID() const
 136 {
 137 //
 138 // Return reference index, using the "GetID" method
 139 // of the possible source object.
 140 //
 141
 142   AliESDtrack *esd = dynamic_cast<AliESDtrack*>(fRef);
 143   if (esd) return esd->GetID();
 144
 145   AliAODTrack *aod = dynamic_cast<AliAODTrack*>(fRef);
 146   if (aod) return aod->GetID();
 147
 148   return GetLabel();
 149 }
 150
 151 //_____________________________________________________________________________
 152 AliPID::EParticleType AliRsnDaughter::RealisticPID() const
 153 {
 154 //
 155 // Return the "realistic" PID of this track,
 156 // i.e. the particle species to which corresponds the largest PID probability.
 157 //
 158
 159   AliPID::EParticleType pid = AliPID::kElectron;
 160   Double_t prob = fPID[0];
 161
 162   Int_t i;
 163   for (i = 1; i < AliPID::kSPECIES; i++) {
 164     if (fPID[i] > prob) {
 165       prob = fPID[i];
 166       pid = (AliPID::EParticleType)i;
 167     }
 168   }
 169
 170   return pid;
 171 }
 172
 173 //_____________________________________________________________________________
 174 AliPID::EParticleType AliRsnDaughter::PerfectPID() const
 175 {
 176 //
 177 // Return the "perfect" PID of this track,
 178 // reading it from the MC information, if available.
 179 //
 180
 181   if (!fParticle) return AliPID::kUnknown;
 182
 183   Int_t absPDG = TMath::Abs(fParticle->GetPdgCode());
 184   switch (absPDG) {
 185     case   11:
 186       return AliPID::kElectron;
 187     case   13:
 188       return AliPID::kMuon;
 189     case  211:
 190       return AliPID::kPion;
 191     case  321:
 192       return AliPID::kKaon;
 193     case 2212:
 194       return AliPID::kProton;
 195     default:
 196       AliDebug(2, Form("PDG code = %d not recognized. Return 'AliPID::kUnknown'", absPDG));
 197       return AliPID::kUnknown;
 198   }
 199 }
 200
 201 //_____________________________________________________________________________
 202 AliPID::EParticleType AliRsnDaughter::PIDType(Double_t &prob) const
 203 {
 204 //
 205 // Return the PID type according to the selected method
 206 // in the argument passed by reference, the probability is stored.
 207 // It will be realistic for realistic PID and 1 for perfect PID.
 208 //
 209
 210   AliPID::EParticleType pid = AssignedPID();
 211
 212   prob = 1.0;
 213   if (fgPIDMethod == kRealistic) prob = fPID[(Int_t)pid];
 214
 215   return pid;
 216 }
 217
 218 //_____________________________________________________________________________
 219 AliPID::EParticleType AliRsnDaughter::AssignedPID() const
 220 {
 221 //
 222 // Return the PID type according to the selected method
 223 // in the argument passed by reference, the probability is stored.
 224 // It will be realistic for realistic PID and 1 for perfect PID.
 225 //
 226
 227   switch (fgPIDMethod) {
 228     case kNoPID:
 229       return AliPID::kUnknown;
 230     case kPerfect:
 231       return PerfectPID();
 232     case kRealistic:
 233       return RealisticPID();
 234     default:
 235       AliWarning("PID method not properly set. Returning realistic PID");
 236       return RealisticPID();
 237   }
 238 }
 239
 240 //_____________________________________________________________________________
 241 Bool_t AliRsnDaughter::CombineWithPriors(Double_t *priors)
 242 {
 243 //
 244 // Combine current PID weights (assumed to be them) with prior probs
 245 //
 246
 247   Int_t       i;
 248   Double_t    sum = 0.0;
 249
 250   // multiply weights and priors
 251   for (i = 0; i < AliPID::kSPECIES; i++) {
 252     fPID[i] = priors[i] * fRef->PID()[i];
 253     sum += fPID[i];
 254   }
 255   if (sum <= (Double_t) 0.) {
 256     AliError(Form("Sum of weights = %f <= 0", sum));
 257     return kFALSE;
 258   }
 259
 260   // normalize
 261   for (i = 0; i < AliPID::kSPECIES; i++) fPID[i] /= sum;
 262
 263   return kTRUE;
 264 }
 265
 266 //_____________________________________________________________________________
 267 AliESDtrack* AliRsnDaughter::GetRefESD()
 268 {
 269 //
 270 // Return a reference in format of ESD track
 271 //
 272
 273   return dynamic_cast<AliESDtrack *>(fRef);
 274 }
 275
 276 //_____________________________________________________________________________
 277 void AliRsnDaughter::FindMotherPDG(AliStack *stack)
 278 {
 279 //
 280 // Searches the stack to find the mother and retrieve its PDG code.
 281 //
 282
 283   if (!stack || !fParticle) return;
 284
 285   Int_t mLabel = fParticle->GetFirstMother();
 286   if (mLabel < 0) {
 287     fMotherPDG = 0;
 288   }
 289   else {
 290     TParticle *mum = stack->Particle(mLabel);
 291     if (mum) fMotherPDG = mum->GetPdgCode();
 292     else fMotherPDG = 0;
 293   }
 294 }
 295
 296 //_____________________________________________________________________________
 297 inline Double_t AliRsnDaughter::GetMCEnergy(Double_t mass)
 298 {
 299 //
 300 // Uses the argument to compute 4-momentum energy
 301 //
 302
 303   if (!fParticle) return 0.0;
 304
 305   Double_t p2 = fParticle->Px()*fParticle->Px();
 306   p2 += fParticle->Py()*fParticle->Py();
 307   p2 += fParticle->Pz()*fParticle->Pz();
 308
 309   return TMath::Sqrt(mass*mass + p2);
 310 }
 311
 312 //_____________________________________________________________________________
 313 void AliRsnDaughter::FindKinkIndex(AliESDtrack *esdTrack)
 314 {
 315 //
 316 // Assign kink index from an ESD track
 317 //
 318
 319   Int_t i, ik[3];
 320   for (i = 0; i < 3; i++) ik[i] = esdTrack->GetKinkIndex(i);
 321
 322   if (ik[0] < 0 || ik[1] < 0 || ik[2] < 0) {
 323     SetKinkMother();
 324   }
 325   else if (ik[0] > 0 || ik[1] > 0 || ik[2] > 0) {
 326     SetKinkDaughter();
 327   }
 328   else SetNoKink();
 329 }
 330
 331 //_____________________________________________________________________________
 332 void AliRsnDaughter::FindKinkIndex(AliAODEvent *event)
 333 {
 334 //
 335 // Assign kink index from an AOD event
 336 //
 337
 338   Int_t iv, id, nD, nV = event->GetNumberOfVertices();
 339   for (iv = 0; iv < nV; iv++) {
 340     AliAODVertex *v = event->GetVertex(iv);
 341     AliAODVertex::AODVtx_t type = (AliAODVertex::AODVtx_t)v->GetType();
 342     if (type != AliAODVertex::kKink) continue;
 343     AliAODTrack *mother = (AliAODTrack*)v->GetParent();
 344     if (mother == (AliAODTrack*)fRef) {
 345       SetKinkMother();
 346       return;
 347     } else {
 348       nD = v->GetNDaughters();
 349       for (id = 0; id < nD; id++) {
 350         AliAODTrack *son = (AliAODTrack*)v->GetDaughter(id);
 351         if (son == (AliAODTrack*)fRef) {
 352           SetKinkDaughter();
 353           return;
 354         }
 355       }
 356     }
 357   }
 358
 359   SetNoKink();
 360 }
 361
 362 //_____________________________________________________________________________
 363 void AliRsnDaughter::Reset()
 364 {
 365 //
 366 // Reset this track to meaningless values
 367 //
 368
 369   fOK = kFALSE;
 370   fKinkIndex = 0;
 371   fParticle = 0x0;
 372   fMotherPDG = 0;
 373   fStatus = 0;
 374   fRef = 0x0;
 375 }
 376
 377 //_____________________________________________________________________________
 378 void AliRsnDaughter::Print(Option_t *option) const
 379 {
 380 //
 381 // Prints the values of data members, using the options:
 382 // - P --> momentum
 383 // - V --> DCA vertex
 384 // - C --> electric charge
 385 // - F --> flags
 386 // - I --> identification (PID, probability and mass)
 387 // - W --> PID weights
 388 // - M --> Montecarlo (from AliRsnMCInfo)
 389 // - L --> index & label
 390 // - A --> angles
 391 // - ALL --> All oprions switched on
 392 //
 393 // Index and label are printed by default.
 394 //
 395
 396   TString opt(option);
 397   opt.ToUpper();
 398
 399   if (opt.Contains("L") || opt.Contains("ALL")) {
 400     cout << ".......Index            : " << GetID() << endl;
 401     cout << ".......Label            : " << GetLabel() << endl;
 402   }
 403   if (opt.Contains("P") || opt.Contains("ALL")) {
 404     cout << ".......Px, Py, Pz, Pt   : " << Px() << ' ' << Py() << ' ' << Pz() << ' ' << Pt() << endl;
 405   }
 406   if (opt.Contains("A") || opt.Contains("ALL")) {
 407     cout << ".......Phi, Theta       : " << Phi() << ' ' << Theta() << endl;
 408   }
 409   if (opt.Contains("V") || opt.Contains("ALL")) {
 410     cout << ".......Vx, Vy, Vz       : " << Xv() << ' ' << Yv() << ' ' << Zv() << endl;
 411   }
 412   if (opt.Contains("I") || opt.Contains("ALL")) {
 413     AliPID::EParticleType type;
 414     Double_t prob;
 415     type = PIDType(prob);
 416     cout << ".......PID & prob       : " << AliPID::ParticleName(type) << ' ' << prob << endl;
 417   }
 418   if (opt.Contains("C") || opt.Contains("ALL")) {
 419     cout << ".......Charge           : " << Charge() << endl;
 420   }
 421   if (opt.Contains("F") || opt.Contains("ALL")) {
 422     cout << ".......Flags            : " << fStatus << endl;
 423   }
 424   if (opt.Contains("W") || opt.Contains("ALL")) {
 425     cout << ".......Weights          : ";
 426     Int_t i;
 427     for (i = 0; i < AliPID::kSPECIES; i++) cout << fPID[i] << ' ';
 428     cout << endl;
 429   }
 430   if (opt.Contains("M") || opt.Contains("ALL")) {
 431     if (fParticle) {
 432       cout << ".......PDG code         : " << fParticle->GetPdgCode() << endl;
 433       cout << ".......Mother (label)   : " << fParticle->GetFirstMother() << endl;
 434       cout << ".......Mother (PDG code): " << fMotherPDG << endl;
 435     } else {
 436       cout << ".......MC info not present" << endl;
 437     }
 438   }
 439 }
 440
 441 //_____________________________________________________________________________
 442 AliPID::EParticleType AliRsnDaughter::InternalType(Int_t pdg)
 443 //
 444 // Return the internal enum value corresponding to the PDG
 445 // code passed as argument, if possible.
 446 // Otherwise, returns 'AliPID::kSPECIES' by default.
 447 //
 448 {
 449   AliPID::EParticleType value;
 450   Int_t absPDG = TMath::Abs(pdg);
 451
 452   switch (absPDG) {
 453     case 11:
 454       value = AliPID::kElectron;
 455       break;
 456     case 13:
 457       value = AliPID::kMuon;
 458       break;
 459     case 211:
 460       value = AliPID::kPion;
 461       break;
 462     case 321:
 463       value = AliPID::kKaon;
 464       break;
 465     case 2212:
 466       value = AliPID::kProton;
 467       break;
 468     default:
 469       value = AliPID::kUnknown;
 470   }
 471   return value;
 472 }