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