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