PWGLF/RESONANCES/AliRsnEvent.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //  This class works as generic interface to an event.
  19 //  Its main purpose is to provide a unique reference which includes all the
  20 //  facilities available in the AliVEvent generic base class, plus all info
  21 //  which could be needed during analysis, which are not in AliVEvent but
  22 //  need to be accessed from ESD or AOD objects, usually in different ways.
  23 //  When MC is available, it is properly taken into account.
  24 //
  25 //  authors: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
  26 //           M. Vala (martin.vala@cern.ch)
  27 //
  28 ////////////////////////////////////////////////////////////////////////////////
  29
  30 #include <Riostream.h>
  31 #include <TArrayF.h>
  32
  33 #include "AliGenEventHeader.h"
  34
  35 #include "AliRsnCutSet.h"
  36 #include "AliRsnEvent.h"
  37
  38 ClassImp(AliRsnEvent)
  39
  40 //_____________________________________________________________________________
  41 AliRsnEvent::AliRsnEvent(AliVEvent *ref, AliVEvent *refMC) :
  42    fRef(ref),
  43    fRefMC(refMC),
  44    fLeading(-1),
  45    fPID(0x0),
  46    fAODList(0x0)
  47 {
  48 //
  49 // Default constructor.
  50 //
  51 }
  52
  53 //_____________________________________________________________________________
  54 AliRsnEvent::AliRsnEvent(const AliRsnEvent &event) :
  55    TObject(event),
  56    fRef(event.fRef),
  57    fRefMC(event.fRefMC),
  58    fLeading(event.fLeading),
  59    fPID(event.fPID),
  60    fAODList(event.fAODList)
  61 {
  62 //
  63 // Copy constructor.
  64 //
  65 }
  66
  67 //_____________________________________________________________________________
  68 AliRsnEvent &AliRsnEvent::operator= (const AliRsnEvent &event)
  69 {
  70 //
  71 // Works in the same way as the copy constructor.
  72 //
  73
  74    TObject::operator=(event);
  75    if (this == &event)
  76       return *this;
  77    fRef             = event.fRef;
  78    fRefMC           = event.fRefMC;
  79    fLeading         = event.fLeading;
  80    fPID             = event.fPID;
  81    fAODList         = event.fAODList;
  82
  83    return (*this);
  84 }
  85
  86 //_____________________________________________________________________________
  87 AliRsnEvent::~AliRsnEvent()
  88 {
  89 //
  90 // Destructor (does nothing since there are not owned pointers)
  91 //
  92 }
  93
  94 //_____________________________________________________________________________
  95 void AliRsnEvent::SetDaughter(AliRsnDaughter &out, Int_t index, Bool_t fromMC)
  96 {
  97 //
  98 // Assigns to the first argument the reference to the i-th track in the ref event.
  99 // What assignment method to be used will depend on the index and on the type of input.
 100 // If the last argument is kTRUE and an MC is referenced, then both fRef and fRefMC will
 101 // point to the MC particle (pure MC analysis)
 102 //
 103
 104    // by default the daughter is reset
 105    // and assigned the used index
 106    out.Reset();
 107    out.SetRsnID(index);
 108    out.SetOwnerEvent(this);
 109
 110    // check input type
 111    if (!InputOK()) return;
 112    Bool_t inputESD = IsESD();
 113
 114    // if it is pure MC, the index tells what particle
 115    // to be read in the stack of MC particles, otherwise
 116    // it is converted into a real collection index
 117    if (fromMC) {
 118       out.SetLabel(index);
 119       Bool_t ok = (inputESD ? SetMCInfoESD(out) : SetMCInfoAOD(out));
 120       if (ok) {
 121          out.SetGood();
 122          out.SetRef(out.GetRefMC());
 123       }
 124    } else {
 125       Int_t trueIndex;
 126       AliRsnDaughter::ERefType type;
 127       if (!ConvertAbsoluteIndex(index, trueIndex, type)) {
 128          AliError(Form("Failed to convert absolute index %d", index));
 129          return;
 130       }
 131       switch (type) {
 132          case AliRsnDaughter::kTrack:
 133             if (inputESD) SetDaughterESDtrack(out, trueIndex); else SetDaughterAODtrack(out, trueIndex);
 134             break;
 135          case AliRsnDaughter::kV0:
 136             if (inputESD) SetDaughterESDv0(out, trueIndex); else SetDaughterAODv0(out, trueIndex);
 137             break;
 138          case AliRsnDaughter::kCascade:
 139             if (inputESD) SetDaughterESDcascade(out, trueIndex); else SetDaughterAODcascade(out, trueIndex);
 140             break;
 141          default:
 142             AliError("Unrecognized daughter type");
 143             return;
 144       }
 145    }
 146 }
 147
 148 //_____________________________________________________________________________
 149 AliRsnDaughter AliRsnEvent::GetDaughter(Int_t i, Bool_t fromMC)
 150 {
 151 //
 152 // Returns a daughter set using same criteria as SetDaughter
 153 //
 154
 155    AliRsnDaughter d;
 156    SetDaughter(d, i, fromMC);
 157    return d;
 158 }
 159
 160 //_____________________________________________________________________________
 161 void AliRsnEvent::SetDaughterESDtrack(AliRsnDaughter &out, Int_t i)
 162 {
 163 //
 164 // Setup the first argument to the track identified by the index.
 165 // When available, adds the MC information and references.
 166 // ---
 167 // Version #1: ESD tracks
 168 //
 169
 170    AliESDEvent *esd = (AliESDEvent *)fRef;
 171
 172    if (i >= 0 && i < esd->GetNumberOfTracks()) {
 173       AliESDtrack *track = esd->GetTrack(i);
 174       if (track) {
 175          out.SetRef(track);
 176          out.SetGood();
 177          // if MC is present, assign label and retrieve corresponding particle
 178          if (fRefMC) {
 179             out.SetLabel(TMath::Abs(track->GetLabel()));
 180             if (!SetMCInfoESD(out)) {
 181                AliWarning("Failed assignment of MC info");
 182             }
 183          }
 184       } else {
 185          AliWarning("Null track");
 186       }
 187    } else {
 188       AliWarning(Form("Overflow: required index = %d, max = %d", i, esd->GetNumberOfTracks()));
 189    }
 190 }
 191
 192 //_____________________________________________________________________________
 193 void AliRsnEvent::SetDaughterAODtrack(AliRsnDaughter &out, Int_t i)
 194 {
 195 //
 196 // Setup the first argument to the track identified by the index.
 197 // When available, adds the MC information and references.
 198 // ---
 199 // Version #2: AOD tracks
 200 //
 201
 202    AliAODEvent *aod = (AliAODEvent *)fRef;
 203
 204    if (i >= 0 && i < aod->GetNumberOfTracks()) {
 205       AliAODTrack *track = dynamic_cast<AliAODTrack*>(aod->GetTrack(i));
 206       if(!track) AliFatal("Not a standard AOD");
 207       if (track) {
 208          out.SetRef(track);
 209          out.SetGood();
 210          // if MC is present, assign label and retrieve corresponding particle
 211          if (fRefMC) {
 212             out.SetLabel(TMath::Abs(track->GetLabel()));
 213             if (!SetMCInfoAOD(out)) {
 214                AliWarning("Failed assignment of MC info");
 215             }
 216          }
 217       } else {
 218          AliWarning("Null track");
 219       }
 220    } else {
 221       AliWarning(Form("Overflow: required index = %d, max = %d", i, aod->GetNumberOfTracks()));
 222    }
 223 }
 224
 225 //_____________________________________________________________________________
 226 void AliRsnEvent::SetDaughterESDv0(AliRsnDaughter &out, Int_t i)
 227 {
 228 //
 229 // Setup the first argument to the track identified by the index.
 230 // When available, adds the MC information and references.
 231 // ---
 232 // Version #3: ESD v0
 233 //
 234
 235    if (i >= 0 && i < fRef->GetNumberOfV0s()) {
 236       AliESDEvent *ev = GetRefESD();
 237       AliESDv0    *v0 = ev->GetV0(i);
 238       if (v0) {
 239          out.SetRef(v0);
 240          out.SetGood();
 241          // if MC is present, retrieve the label of V0 from those of daughters
 242          if (fRefMC) {
 243             AliMCEvent  *mc = (AliMCEvent *)fRefMC;
 244             AliESDtrack *tp = ev->GetTrack(v0->GetPindex());
 245             AliESDtrack *tn = ev->GetTrack(v0->GetNindex());
 246             if (mc && tp && tn) {
 247                Int_t lp = TMath::Abs(tp->GetLabel());
 248                Int_t ln = TMath::Abs(tn->GetLabel());
 249                TParticle *pp = mc->Stack()->Particle(lp);
 250                TParticle *pn = mc->Stack()->Particle(ln);
 251                if (pp && pn) {
 252                   // if their first mothers are the same, the V0 is true
 253                   // otherwise label remains '-1' --> fake V0
 254                   if (pp->GetFirstMother() == pn->GetFirstMother() && pp->GetFirstMother() >= 0) {
 255                      out.SetLabel(pp->GetFirstMother());
 256                      //patch for k0s/k0l
 257                      TParticle *mom = mc->Stack()->Particle(pn->GetFirstMother());
 258                      if(mom->GetPdgCode() == 310) {
 259                         //take the mother of the k0s which is a k0 (311)
 260                         out.SetLabel(mom->GetFirstMother());
 261                      }
 262                      SetMCInfoESD(out);
 263                   }
 264                }
 265             }
 266          }
 267       }
 268    }
 269 }
 270
 271 //_____________________________________________________________________________
 272 void AliRsnEvent::SetDaughterAODv0(AliRsnDaughter &out, Int_t i)
 273 {
 274 //
 275 // Setup the first argument to the track identified by the index.
 276 // When available, adds the MC information and references.
 277 // ---
 278 // Version #4: AOD v0
 279 //
 280
 281    if (i >= 0 && i < fRef->GetNumberOfV0s()) {
 282       AliAODEvent *ev = (AliAODEvent *)fRef;
 283       AliAODv0    *v0 = ev->GetV0(i);
 284       if (v0) {
 285          out.SetRef(v0);
 286          out.SetGood();
 287          if (fRefMC) {
 288             AliAODEvent  *mc = (AliAODEvent *)fRefMC;
 289             TClonesArray *mcArray = (TClonesArray *)mc->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 290             AliAODTrack  *tp  = (AliAODTrack *)v0->GetDaughter(0);
 291             AliAODTrack  *tn  = (AliAODTrack *)v0->GetDaughter(1);
 292             if (mcArray && tp && tn) {
 293                Int_t lp = TMath::Abs(tp->GetLabel());
 294                Int_t ln = TMath::Abs(tn->GetLabel());
 295                AliAODMCParticle *pp = (AliAODMCParticle *)mcArray->At(lp);
 296                AliAODMCParticle *pn = (AliAODMCParticle *)mcArray->At(ln);
 297                if (pp && pn) {
 298                   // if their first mothers are the same, the V0 is true
 299                   // otherwise label remains '-1' --> fake V0
 300                   if (pp->GetMother() == pn->GetMother() && pp->GetMother() >= 0) {
 301                      out.SetLabel(pp->GetMother());
 302                      SetMCInfoAOD(out);
 303                   }
 304                }
 305             }
 306          }
 307       }
 308    }
 309 }
 310
 311 //_____________________________________________________________________________
 312 void AliRsnEvent::SetDaughterESDcascade(AliRsnDaughter &out, Int_t i)
 313 {
 314 //
 315 // Setup the first argument to the track identified by the index.
 316 // When available, adds the MC information and references.
 317 // ---
 318 // Version #3: ESD cascade
 319 //
 320
 321    if (i >= 0 && i < fRef->GetNumberOfCascades()) {
 322       AliESDEvent   *ev   = GetRefESD();
 323       AliESDcascade *casc = ev->GetCascade(i);
 324       if (casc) {
 325          out.SetRef(casc);
 326          out.SetGood();
 327          if (fRefMC) {
 328
 329          }
 330       }
 331    }
 332 }
 333
 334 //_____________________________________________________________________________
 335 void AliRsnEvent::SetDaughterAODcascade(AliRsnDaughter &out, Int_t i)
 336 {
 337 //
 338 // Setup the first argument to the track identified by the index.
 339 // When available, adds the MC information and references.
 340 // ---
 341 // Version #4: AOD cascade
 342 //
 343
 344    if (i >= 0 && i < fRef->GetNumberOfCascades()) {
 345       AliAODEvent *ev = GetRefAOD();
 346       AliAODv0    *casc = ev->GetCascade(i);
 347       if (casc) {
 348          out.SetRef(casc);
 349          out.SetGood();
 350          if (fRefMC) {
 351
 352          }
 353       }
 354    }
 355 }
 356
 357 //_____________________________________________________________________________
 358 Bool_t AliRsnEvent::SetMCInfoESD(AliRsnDaughter &out)
 359 {
 360 //
 361 // Using the label assigned to the daughter, searches for the MC informations:
 362 // -- MC reference
 363 // -- mother
 364 //
 365
 366    // if label makes no sense --> failed
 367    Int_t label = out.GetLabel();
 368    if (label < 0 || !fRefMC) return kFALSE;
 369
 370    // get number of particles
 371    Int_t nMC = fRefMC->GetNumberOfTracks();
 372
 373    // if label too large --> failed
 374    if (label >= nMC) {
 375      AliDebug(4, Form("Stack overflow: track label = %d -- stack maximum = %d", label, nMC));
 376       return kFALSE;
 377    }
 378
 379    // retrieve particle
 380    AliMCEvent    *mc = (AliMCEvent *)fRefMC;
 381    AliMCParticle *mcPart = (AliMCParticle *)mc->GetTrack(label);
 382
 383    // if particle = NULL --> failed
 384    if (!mcPart) {
 385       AliWarning(Form("Stack discontinuity: label %d refers to a NULL object", label));
 386       return kFALSE;
 387    }
 388    // otherwise --> success
 389    out.SetRefMC(mcPart);
 390
 391    // if the particle is not primary, find the mother and get its PDG
 392    Int_t imum = mcPart->Particle()->GetFirstMother();
 393    if (imum >= 0 && imum < nMC) {
 394       AliMCParticle *mcMother = (AliMCParticle *)mc->GetTrack(imum);
 395       if (mcMother) {
 396          out.SetMotherPDG(mcMother->Particle()->GetPdgCode());
 397       } else {
 398          AliWarning(Form("Stack discontinuity: label mother %d refers to a NULL object", imum));
 399       }
 400    } else {
 401      AliDebug(4, Form("Stack overflow: mother label = %d -- stack maximum = %d", imum, nMC));
 402    }
 403
 404    return kTRUE;
 405 }
 406
 407 //_____________________________________________________________________________
 408 Bool_t AliRsnEvent::SetMCInfoAOD(AliRsnDaughter &out)
 409 {
 410 //
 411 // Using the label assigned to the daughter, searches for the MC informations:
 412 // -- MC reference
 413 // -- mother
 414 //
 415
 416    // if label makes no sense --> failed
 417    Int_t label = out.GetLabel();
 418    if (label < 0 || !fRefMC) return kFALSE;
 419
 420    // retrieve particle
 421    AliAODEvent  *mc = (AliAODEvent *)fRefMC;
 422    TClonesArray *mcArray = (TClonesArray *)mc->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 423
 424    // get number of particles
 425    Int_t nMC = mcArray->GetEntriesFast();
 426
 427    // if label too large --> failed
 428    if (label >= nMC) {
 429      AliDebug(4, Form("Stack overflow: track label = %d -- stack maximum = %d", label, nMC));
 430       return kFALSE;
 431    }
 432
 433    // if particle = NULL --> failed
 434    AliAODMCParticle *mcPart = (AliAODMCParticle *)mcArray->At(label);
 435    if (!mcPart) {
 436       AliWarning(Form("Stack discontinuity: label %d refers to a NULL object", label));
 437       return kFALSE;
 438    }
 439    // otherwise --> success
 440    out.SetRefMC(mcPart);
 441
 442    // if the particle is not primary, find the mother and get its PDG
 443    Int_t imum = mcPart->GetMother();
 444    if (imum >= 0 && imum < nMC) {
 445       AliAODMCParticle *mcMother = (AliAODMCParticle *)mcArray->At(imum);
 446       if (mcMother) {
 447          out.SetMotherPDG(mcMother->GetPdgCode());
 448       } else {
 449          AliWarning(Form("Stack discontinuity: label mother %d refers to a NULL object", imum));
 450       }
 451    } else if (imum >= nMC) {
 452      AliDebug(4, Form("Stack overflow: mother label = %d -- stack maximum = %d", imum, nMC));
 453    }
 454
 455    return kTRUE;
 456 }
 457
 458 //_____________________________________________________________________________
 459 Bool_t AliRsnEvent::ConvertAbsoluteIndex(Int_t index, Int_t &realIndex, AliRsnDaughter::ERefType &type)
 460 {
 461 //
 462 // Using the phylosophy of the absolute index, which loops over
 463 // all tracks, V0s and cascades, returns the result of a check
 464 // on it (first argument) based on this criterion:
 465 // 1) if the absolute ID is smaller than number of tracks,
 466 //    return itself and the type 'track'
 467 // 2) if the absolute ID is larger than number of tracks, subtract it
 468 //    and if the result is smaller than number of V0s,
 469 //    return the corresponding V0 index and type
 470 // 3) if the absolute ID is larger than number of tracks + V0s, subtract them
 471 //    and if the result is smaller than number of cascades,
 472 //    return the corresponding cascade index and type
 473 // The results of this check are stored in the reference arguments, while the outcome of
 474 // the function is kTRUE if one of these checks was successful, otherwise it is kFALSE,
 475 // meaning that the absolute index reached the end.
 476 //
 477
 478    Int_t nTracks   = fRef->GetNumberOfTracks();
 479    Int_t nV0s      = fRef->GetNumberOfV0s();
 480    Int_t nCascades = fRef->GetNumberOfCascades();
 481
 482    if (index < nTracks) {
 483       realIndex = index;
 484       type = AliRsnDaughter::kTrack;
 485       return kTRUE;
 486    } else if (index >= nTracks && index < nTracks + nV0s) {
 487       realIndex = index - nTracks;
 488       type = AliRsnDaughter::kV0;
 489       return kTRUE;
 490    } else if (index >= nTracks + nV0s && index < nTracks + nV0s + nCascades) {
 491       realIndex = index - nTracks - nV0s;
 492       type = AliRsnDaughter::kCascade;
 493       return kTRUE;
 494    }
 495
 496    realIndex = -1;
 497    type = AliRsnDaughter::kNoType;
 498    return kFALSE;
 499 }
 500
 501 //_____________________________________________________________________________
 502 Int_t AliRsnEvent::ConvertRealIndex(Int_t index, AliRsnDaughter::ERefType type)
 503 {
 504 //
 505 // Translates a pair made by index + object type into the corresponding
 506 // absolute index, which is set to -1 in case the real index overflows.
 507 //
 508
 509    Int_t nTracks   = fRef->GetNumberOfTracks();
 510    Int_t nV0s      = fRef->GetNumberOfV0s();
 511    Int_t nCascades = fRef->GetNumberOfCascades();
 512
 513    switch (type) {
 514       case AliRsnDaughter::kTrack:
 515          if (index >= 0 && index < nTracks)
 516             return index;
 517          else
 518             return -1;
 519       case AliRsnDaughter::kV0:
 520          if (index >= 0 && index < nV0s)
 521             return nTracks + index;
 522          else
 523             return -1;
 524       case AliRsnDaughter::kCascade:
 525          if (index >= 0 && index < nCascades)
 526             return nTracks + nV0s + index;
 527          else
 528             return -1;
 529       default:
 530          return -1;
 531    }
 532 }
 533
 534 //_____________________________________________________________________________
 535 Int_t AliRsnEvent::SelectLeadingParticle(AliRsnCutSet *cuts)
 536 {
 537 //
 538 // Searches the collection of all particles with given PID type and charge,
 539 // and returns the one with largest momentum, provided that it is greater than 1st argument.
 540 // If one specifies AliRsnPID::kUnknown as type or AliRsnDaughter::kNoPID as method,
 541 // the check is done over all particles irrespectively of their PID.
 542 // If the sign argument is '+' or '-', the check is done over the particles of that charge,
 543 // otherwise it is done irrespectively of the charge.
 544 //
 545
 546    // check input type
 547    Bool_t inputESD = IsESD();
 548    if (!inputESD && !IsAOD()) {
 549       AliError("Need to process ESD or AOD input");
 550       return -1;
 551    }
 552
 553    Double_t ptMax = 0.0;
 554    Int_t i, nTracks = fRef->GetNumberOfTracks();
 555
 556    fLeading = -1;
 557    AliRsnDaughter leading;
 558
 559    for (i = 0; i < nTracks; i++) {
 560       if (inputESD)
 561          SetDaughterESDtrack(leading, i);
 562       else
 563          SetDaughterAODtrack(leading, i);
 564       if (!leading.IsOK()) {
 565          AliDebugClass(1, Form("Failed assignment of track %d", i));
 566          continue;
 567       }
 568       if (cuts && !cuts->IsSelected(&leading)) {
 569          AliDebugClass(1, Form("Track %d didn't pass cuts", i));
 570          continue;
 571       }
 572       // check if it has largest momentum
 573       if (leading.GetRef()->Pt() > ptMax) {
 574          ptMax = leading.GetRef()->Pt();
 575          fLeading = i;
 576       }
 577    }
 578
 579    return fLeading;
 580 }