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 = aod->GetTrack(i);
 206       if (track) {
 207          out.SetRef(track);
 208          out.SetGood();
 209          // if MC is present, assign label and retrieve corresponding particle
 210          if (fRefMC) {
 211             out.SetLabel(TMath::Abs(track->GetLabel()));
 212             if (!SetMCInfoAOD(out)) {
 213                AliWarning("Failed assignment of MC info");
 214             }
 215          }
 216       } else {
 217          AliWarning("Null track");
 218       }
 219    } else {
 220       AliWarning(Form("Overflow: required index = %d, max = %d", i, aod->GetNumberOfTracks()));
 221    }
 222 }
 223
 224 //_____________________________________________________________________________
 225 void AliRsnEvent::SetDaughterESDv0(AliRsnDaughter &out, Int_t i)
 226 {
 227 //
 228 // Setup the first argument to the track identified by the index.
 229 // When available, adds the MC information and references.
 230 // ---
 231 // Version #3: ESD v0
 232 //
 233
 234    if (i >= 0 && i < fRef->GetNumberOfV0s()) {
 235       AliESDEvent *ev = GetRefESD();
 236       AliESDv0    *v0 = ev->GetV0(i);
 237       if (v0) {
 238          out.SetRef(v0);
 239          out.SetGood();
 240          // if MC is present, retrieve the label of V0 from those of daughters
 241          if (fRefMC) {
 242             AliMCEvent  *mc = (AliMCEvent *)fRefMC;
 243             AliESDtrack *tp = ev->GetTrack(v0->GetPindex());
 244             AliESDtrack *tn = ev->GetTrack(v0->GetNindex());
 245             if (mc && tp && tn) {
 246                Int_t lp = TMath::Abs(tp->GetLabel());
 247                Int_t ln = TMath::Abs(tn->GetLabel());
 248                TParticle *pp = mc->Stack()->Particle(lp);
 249                TParticle *pn = mc->Stack()->Particle(ln);
 250                if (pp && pn) {
 251                   // if their first mothers are the same, the V0 is true
 252                   // otherwise label remains '-1' --> fake V0
 253                   if (pp->GetFirstMother() == pn->GetFirstMother() && pp->GetFirstMother() >= 0) {
 254                      out.SetLabel(pp->GetFirstMother());
 255                      SetMCInfoESD(out);
 256                   }
 257                }
 258             }
 259          }
 260       }
 261    }
 262 }
 263
 264 //_____________________________________________________________________________
 265 void AliRsnEvent::SetDaughterAODv0(AliRsnDaughter &out, Int_t i)
 266 {
 267 //
 268 // Setup the first argument to the track identified by the index.
 269 // When available, adds the MC information and references.
 270 // ---
 271 // Version #4: AOD v0
 272 //
 273
 274    if (i >= 0 && i < fRef->GetNumberOfV0s()) {
 275       AliAODEvent *ev = (AliAODEvent *)fRef;
 276       AliAODv0    *v0 = ev->GetV0(i);
 277       if (v0) {
 278          out.SetRef(v0);
 279          out.SetGood();
 280          if (fRefMC) {
 281             AliAODEvent  *mc = (AliAODEvent *)fRefMC;
 282             TClonesArray *mcArray = (TClonesArray *)mc->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 283             AliAODTrack  *tp  = (AliAODTrack *)v0->GetDaughter(0);
 284             AliAODTrack  *tn  = (AliAODTrack *)v0->GetDaughter(1);
 285             if (mcArray && tp && tn) {
 286                Int_t lp = TMath::Abs(tp->GetLabel());
 287                Int_t ln = TMath::Abs(tn->GetLabel());
 288                AliAODMCParticle *pp = (AliAODMCParticle *)mcArray->At(lp);
 289                AliAODMCParticle *pn = (AliAODMCParticle *)mcArray->At(ln);
 290                if (pp && pn) {
 291                   // if their first mothers are the same, the V0 is true
 292                   // otherwise label remains '-1' --> fake V0
 293                   if (pp->GetMother() == pn->GetMother() && pp->GetMother() >= 0) {
 294                      out.SetLabel(pp->GetMother());
 295                      SetMCInfoAOD(out);
 296                   }
 297                }
 298             }
 299          }
 300       }
 301    }
 302 }
 303
 304 //_____________________________________________________________________________
 305 void AliRsnEvent::SetDaughterESDcascade(AliRsnDaughter &out, Int_t i)
 306 {
 307 //
 308 // Setup the first argument to the track identified by the index.
 309 // When available, adds the MC information and references.
 310 // ---
 311 // Version #3: ESD cascade
 312 //
 313
 314    if (i >= 0 && i < fRef->GetNumberOfCascades()) {
 315       AliESDEvent   *ev   = GetRefESD();
 316       AliESDcascade *casc = ev->GetCascade(i);
 317       if (casc) {
 318          out.SetRef(casc);
 319          out.SetGood();
 320          if (fRefMC) {
 321
 322          }
 323       }
 324    }
 325 }
 326
 327 //_____________________________________________________________________________
 328 void AliRsnEvent::SetDaughterAODcascade(AliRsnDaughter &out, Int_t i)
 329 {
 330 //
 331 // Setup the first argument to the track identified by the index.
 332 // When available, adds the MC information and references.
 333 // ---
 334 // Version #4: AOD cascade
 335 //
 336
 337    if (i >= 0 && i < fRef->GetNumberOfCascades()) {
 338       AliAODEvent *ev = GetRefAOD();
 339       AliAODv0    *casc = ev->GetCascade(i);
 340       if (casc) {
 341          out.SetRef(casc);
 342          out.SetGood();
 343          if (fRefMC) {
 344
 345          }
 346       }
 347    }
 348 }
 349
 350 //_____________________________________________________________________________
 351 Bool_t AliRsnEvent::SetMCInfoESD(AliRsnDaughter &out)
 352 {
 353 //
 354 // Using the label assigned to the daughter, searches for the MC informations:
 355 // -- MC reference
 356 // -- mother
 357 //
 358
 359    // if label makes no sense --> failed
 360    Int_t label = out.GetLabel();
 361    if (label < 0 || !fRefMC) return kFALSE;
 362
 363    // get number of particles
 364    Int_t nMC = fRefMC->GetNumberOfTracks();
 365
 366    // if label too large --> failed
 367    if (label >= nMC) {
 368       AliWarning(Form("Stack overflow: track label = %d -- stack maximum = %d", label, nMC));
 369       return kFALSE;
 370    }
 371
 372    // retrieve particle
 373    AliMCEvent    *mc = (AliMCEvent *)fRefMC;
 374    AliMCParticle *mcPart = (AliMCParticle *)mc->GetTrack(label);
 375
 376    // if particle = NULL --> failed
 377    if (!mcPart) {
 378       AliWarning(Form("Stack discontinuity: label %d refers to a NULL object", label));
 379       return kFALSE;
 380    }
 381    // otherwise --> success
 382    out.SetRefMC(mcPart);
 383
 384    // if the particle is not primary, find the mother and get its PDG
 385    Int_t imum = mcPart->Particle()->GetFirstMother();
 386    if (imum >= 0 && imum < nMC) {
 387       AliMCParticle *mcMother = (AliMCParticle *)mc->GetTrack(imum);
 388       if (mcMother) {
 389          out.SetMotherPDG(mcMother->Particle()->GetPdgCode());
 390       } else {
 391          AliWarning(Form("Stack discontinuity: label mother %d refers to a NULL object", imum));
 392       }
 393    } else {
 394       AliWarning(Form("Stack overflow: mother label = %d -- stack maximum = %d", imum, nMC));
 395    }
 396
 397    return kTRUE;
 398 }
 399
 400 //_____________________________________________________________________________
 401 Bool_t AliRsnEvent::SetMCInfoAOD(AliRsnDaughter &out)
 402 {
 403 //
 404 // Using the label assigned to the daughter, searches for the MC informations:
 405 // -- MC reference
 406 // -- mother
 407 //
 408
 409    // if label makes no sense --> failed
 410    Int_t label = out.GetLabel();
 411    if (label < 0 || !fRefMC) return kFALSE;
 412
 413    // retrieve particle
 414    AliAODEvent  *mc = (AliAODEvent *)fRefMC;
 415    TClonesArray *mcArray = (TClonesArray *)mc->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 416
 417    // get number of particles
 418    Int_t nMC = mcArray->GetEntriesFast();
 419
 420    // if label too large --> failed
 421    if (label >= nMC) {
 422       AliWarning(Form("Stack overflow: track label = %d -- stack maximum = %d", label, nMC));
 423       return kFALSE;
 424    }
 425
 426    // if particle = NULL --> failed
 427    AliAODMCParticle *mcPart = (AliAODMCParticle *)mcArray->At(label);
 428    if (!mcPart) {
 429       AliWarning(Form("Stack discontinuity: label %d refers to a NULL object", label));
 430       return kFALSE;
 431    }
 432    // otherwise --> success
 433    out.SetRefMC(mcPart);
 434
 435    // if the particle is not primary, find the mother and get its PDG
 436    Int_t imum = mcPart->GetMother();
 437    if (imum >= 0 && imum < nMC) {
 438       AliAODMCParticle *mcMother = (AliAODMCParticle *)mcArray->At(imum);
 439       if (mcMother) {
 440          out.SetMotherPDG(mcMother->GetPdgCode());
 441       } else {
 442          AliWarning(Form("Stack discontinuity: label mother %d refers to a NULL object", imum));
 443       }
 444    } else if (imum >= nMC) {
 445       AliWarning(Form("Stack overflow: mother label = %d -- stack maximum = %d", imum, nMC));
 446    }
 447
 448    return kTRUE;
 449 }
 450
 451 //_____________________________________________________________________________
 452 Bool_t AliRsnEvent::ConvertAbsoluteIndex(Int_t index, Int_t &realIndex, AliRsnDaughter::ERefType &type)
 453 {
 454 //
 455 // Using the phylosophy of the absolute index, which loops over
 456 // all tracks, V0s and cascades, returns the result of a check
 457 // on it (first argument) based on this criterion:
 458 // 1) if the absolute ID is smaller than number of tracks,
 459 //    return itself and the type 'track'
 460 // 2) if the absolute ID is larger than number of tracks, subtract it
 461 //    and if the result is smaller than number of V0s,
 462 //    return the corresponding V0 index and type
 463 // 3) if the absolute ID is larger than number of tracks + V0s, subtract them
 464 //    and if the result is smaller than number of cascades,
 465 //    return the corresponding cascade index and type
 466 // The results of this check are stored in the reference arguments, while the outcome of
 467 // the function is kTRUE if one of these checks was successful, otherwise it is kFALSE,
 468 // meaning that the absolute index reached the end.
 469 //
 470
 471    Int_t nTracks   = fRef->GetNumberOfTracks();
 472    Int_t nV0s      = fRef->GetNumberOfV0s();
 473    Int_t nCascades = fRef->GetNumberOfCascades();
 474
 475    if (index < nTracks) {
 476       realIndex = index;
 477       type = AliRsnDaughter::kTrack;
 478       return kTRUE;
 479    } else if (index >= nTracks && index < nTracks + nV0s) {
 480       realIndex = index - nTracks;
 481       type = AliRsnDaughter::kV0;
 482       return kTRUE;
 483    } else if (index >= nTracks + nV0s && index < nTracks + nV0s + nCascades) {
 484       realIndex = index - nTracks - nV0s;
 485       type = AliRsnDaughter::kCascade;
 486       return kTRUE;
 487    }
 488
 489    realIndex = -1;
 490    type = AliRsnDaughter::kNoType;
 491    return kFALSE;
 492 }
 493
 494 //_____________________________________________________________________________
 495 Int_t AliRsnEvent::ConvertRealIndex(Int_t index, AliRsnDaughter::ERefType type)
 496 {
 497 //
 498 // Translates a pair made by index + object type into the corresponding
 499 // absolute index, which is set to -1 in case the real index overflows.
 500 //
 501
 502    Int_t nTracks   = fRef->GetNumberOfTracks();
 503    Int_t nV0s      = fRef->GetNumberOfV0s();
 504    Int_t nCascades = fRef->GetNumberOfCascades();
 505
 506    switch (type) {
 507       case AliRsnDaughter::kTrack:
 508          if (index >= 0 && index < nTracks)
 509             return index;
 510          else
 511             return -1;
 512       case AliRsnDaughter::kV0:
 513          if (index >= 0 && index < nV0s)
 514             return nTracks + index;
 515          else
 516             return -1;
 517       case AliRsnDaughter::kCascade:
 518          if (index >= 0 && index < nCascades)
 519             return nTracks + nV0s + index;
 520          else
 521             return -1;
 522       default:
 523          return -1;
 524    }
 525 }
 526
 527 //_____________________________________________________________________________
 528 Int_t AliRsnEvent::SelectLeadingParticle(AliRsnCutSet *cuts)
 529 {
 530 //
 531 // Searches the collection of all particles with given PID type and charge,
 532 // and returns the one with largest momentum, provided that it is greater than 1st argument.
 533 // If one specifies AliRsnPID::kUnknown as type or AliRsnDaughter::kNoPID as method,
 534 // the check is done over all particles irrespectively of their PID.
 535 // If the sign argument is '+' or '-', the check is done over the particles of that charge,
 536 // otherwise it is done irrespectively of the charge.
 537 //
 538
 539    // check input type
 540    Bool_t inputESD = IsESD();
 541    if (!inputESD && !IsAOD()) {
 542       AliError("Need to process ESD or AOD input");
 543       return -1;
 544    }
 545
 546    Double_t ptMax = 0.0;
 547    Int_t i, nTracks = fRef->GetNumberOfTracks();
 548
 549    fLeading = -1;
 550    AliRsnDaughter leading;
 551
 552    for (i = 0; i < nTracks; i++) {
 553       if (inputESD)
 554          SetDaughterESDtrack(leading, i);
 555       else
 556          SetDaughterAODtrack(leading, i);
 557       if (!leading.IsOK()) {
 558          AliDebugClass(1, Form("Failed assignment of track %d", i));
 559          continue;
 560       }
 561       if (cuts && !cuts->IsSelected(&leading)) {
 562          AliDebugClass(1, Form("Track %d didn't pass cuts", i));
 563          continue;
 564       }
 565       // check if it has largest momentum
 566       if (leading.GetRef()->Pt() > ptMax) {
 567          ptMax = leading.GetRef()->Pt();
 568          fLeading = i;
 569       }
 570    }
 571
 572    return fLeading;
 573 }