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