AliParticleContainer.cxx

Go to the documentation of this file.

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
#include <algorithm>
#include <iostream>
#include <vector>
#include <TClonesArray.h>

#include "AliVEvent.h"
#include "AliLog.h"

#include "AliTLorentzVector.h"
#include "AliParticleContainer.h"

ClassImp(AliParticleContainer)

AliParticleContainer::AliParticleContainer():
  AliEmcalContainer(),
  fMinDistanceTPCSectorEdge(-1),
  fCharge(-1),
  fGeneratorIndex(-1)
{
  fBaseClassName = "AliVParticle";
  SetClassName("AliVParticle");
}

AliParticleContainer::AliParticleContainer(const char *name) :
  AliEmcalContainer(name),
  fMinDistanceTPCSectorEdge(-1),
  fCharge(-1),
  fGeneratorIndex(-1)
{
  fBaseClassName = "AliVParticle";
  SetClassName("AliVParticle");
}

AliVParticle* AliParticleContainer::GetLeadingParticle(const char* opt) 
{
  TString option(opt);
  option.ToLower();

  Int_t tempID = fCurrentID;
  ResetCurrentID();

  AliVParticle *partMax = GetNextAcceptParticle();
  AliVParticle *part = 0;

  if (option.Contains("p")) {
    while ((part = GetNextAcceptParticle())) {
      if (part->P() > partMax->P()) partMax = part;
    }
  }
  else {
    while ((part = GetNextAcceptParticle())) {
      if (part->Pt() > partMax->Pt()) partMax = part;
    }
  }

  fCurrentID = tempID;

  return partMax;
}

AliVParticle* AliParticleContainer::GetParticle(Int_t i) const 
{
  //
  if (i == -1) i = fCurrentID;
  if (i < 0 || i >= this->fClArray->GetEntriesFast()) return 0;
  AliVParticle *vp = static_cast<AliVParticle*>(fClArray->At(i));
  return vp;
}

AliVParticle* AliParticleContainer::GetAcceptParticle(Int_t i) const
{
  UInt_t rejectionReason = 0;
  if (i == -1) i = fCurrentID;
  if (AcceptParticle(i, rejectionReason)) {
      return GetParticle(i);
  }
  else {
    AliDebug(2,"Particle not accepted.");
    return 0;
  }
}

AliVParticle* AliParticleContainer::GetNextAcceptParticle()
{
  const Int_t n = GetNEntries();
  AliVParticle *p = 0;
  do {
    fCurrentID++;
    if (fCurrentID >= n) break;
    p = GetAcceptParticle(fCurrentID);
  } while (!p);

  return p;
}

AliVParticle* AliParticleContainer::GetNextParticle()
{
  const Int_t n = GetNEntries();
  AliVParticle *p = 0;
  do {
    fCurrentID++;
    if (fCurrentID >= n) break;
    p = GetParticle(fCurrentID);
  } while (!p);

  return p;
}

Bool_t AliParticleContainer::GetMomentumFromParticle(TLorentzVector &mom, const AliVParticle* part, Double_t mass) const
{
  if (part) {
    if (mass < 0) mass = part->M();
    mom.SetPtEtaPhiM(part->Pt(), part->Eta(), part->Phi(), mass);
    return kTRUE;
  }
  else {
    mom.SetPtEtaPhiM(0, 0, 0, 0);
    return kFALSE;
  }
}

Bool_t AliParticleContainer::GetMomentumFromParticle(TLorentzVector &mom, const AliVParticle* part) const
{
  return GetMomentumFromParticle(mom,part,fMassHypothesis);
}

Bool_t AliParticleContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
{
  if (i == -1) i = fCurrentID;
  AliVParticle *vp = GetParticle(i);
  return GetMomentumFromParticle(mom, vp);
}

Bool_t AliParticleContainer::GetNextMomentum(TLorentzVector &mom)
{
  AliVParticle *vp = GetNextParticle();
  return GetMomentumFromParticle(mom, vp);
}

Bool_t AliParticleContainer::GetAcceptMomentum(TLorentzVector &mom, Int_t i) const
{
  if (i == -1) i = fCurrentID;
  AliVParticle *vp = GetAcceptParticle(i);
  return GetMomentumFromParticle(mom, vp);
}

Bool_t AliParticleContainer::GetNextAcceptMomentum(TLorentzVector &mom)
{
  AliVParticle *vp = GetNextAcceptParticle();
  return GetMomentumFromParticle(mom, vp);
}

Bool_t AliParticleContainer::AcceptParticle(const AliVParticle *vp, UInt_t &rejectionReason) const
{
  Bool_t r = ApplyParticleCuts(vp, rejectionReason);
  if (!r) return kFALSE;

  AliTLorentzVector mom;

  Int_t id = fClArray->IndexOf(vp);
  if (id >= 0) {
    GetMomentum(mom, id);
  }
  else {
    GetMomentumFromParticle(mom, vp);
  }

  return ApplyKinematicCuts(mom, rejectionReason);
}

Bool_t AliParticleContainer::AcceptParticle(Int_t i, UInt_t &rejectionReason) const
{
  Bool_t r = ApplyParticleCuts(GetParticle(i), rejectionReason);
  if (!r) return kFALSE;

  AliTLorentzVector mom;
  GetMomentum(mom, i);

  return ApplyKinematicCuts(mom, rejectionReason);
}

Bool_t AliParticleContainer::ApplyParticleCuts(const AliVParticle* vp, UInt_t &rejectionReason) const
{
  if (!vp) {
    rejectionReason |= kNullObject;
    return kFALSE;
  }

  if (vp->TestBits(fBitMap) != (Int_t)fBitMap) {
    rejectionReason |= kBitMapCut;
    return kFALSE;
  }

  if (fMinMCLabel >= 0 && TMath::Abs(vp->GetLabel()) > fMinMCLabel) {
    rejectionReason |= kMCLabelCut;
    return kFALSE;
  }

  if (fMaxMCLabel >= 0 && TMath::Abs(vp->GetLabel()) < fMaxMCLabel) {
    rejectionReason |= kMCLabelCut;
    return kFALSE;
  }

  if (fCharge>=0 && fCharge != vp->Charge()) {
    rejectionReason |= kChargeCut;
    return kFALSE;
  }

  if (fGeneratorIndex >= 0 && fGeneratorIndex != vp->GetGeneratorIndex()) {
    rejectionReason |= kMCGeneratorCut;
    return kFALSE;
  }

  return kTRUE;
}

Bool_t AliParticleContainer::ApplyKinematicCuts(const AliTLorentzVector& mom, UInt_t &rejectionReason) const
{
  if(fMinDistanceTPCSectorEdge>0.) {
    const Double_t pi = TMath::Pi();
    const Double_t kSector = pi/9;
    Double_t phiDist = TMath::Abs(mom.Phi() - TMath::FloorNint(mom.Phi()/kSector)*kSector);
    if(phiDist<fMinDistanceTPCSectorEdge) {
      rejectionReason |= kMinDistanceTPCSectorEdgeCut;
      return kFALSE;
    }
  }

  return AliEmcalContainer::ApplyKinematicCuts(mom, rejectionReason);
}

Int_t AliParticleContainer::GetNAcceptedParticles() const
{
  Int_t nPart = 0;
  for(int ipart = 0; ipart < this->GetNParticles(); ipart++){
    UInt_t rejectionReason = 0;
    if(this->AcceptParticle(ipart, rejectionReason)) nPart++;
  }
  return nPart;
}

const char* AliParticleContainer::GetTitle() const
{
  static TString trackString;

  if (GetMinPt() == 0) {
    trackString = TString::Format("%s_pT0000", GetArrayName().Data());
  }
  else if (GetMinPt() < 1.0) {
    trackString = TString::Format("%s_pT0%3.0f", GetArrayName().Data(), GetMinPt()*1000.0);
  }
  else {
    trackString = TString::Format("%s_pT%4.0f", GetArrayName().Data(), GetMinPt()*1000.0);
  }

  return trackString.Data();
}

AliParticleContainer::all_iterator AliParticleContainer::begin() const {
  return all_iterator(this, 0, true);
}

AliParticleContainer::all_iterator AliParticleContainer::end() const {
  return all_iterator(this, this->GetNParticles(), true);
}

AliParticleContainer::all_iterator AliParticleContainer::rbegin() const {
  return all_iterator(this, this->GetNParticles(), false);
}

AliParticleContainer::all_iterator AliParticleContainer::rend() const {
  return all_iterator(this, -1, false);
}
AliParticleContainer::accept_iterator AliParticleContainer::accept_begin() const{
  return accept_iterator(this, 0, true);
}

AliParticleContainer::accept_iterator AliParticleContainer::accept_end() const{
  return accept_iterator(this, this->GetNAcceptedParticles(), true);
}

AliParticleContainer::accept_iterator AliParticleContainer::accept_rbegin() const{
  return accept_iterator(this, this->GetNAcceptedParticles()-1, false);
}

AliParticleContainer::accept_iterator AliParticleContainer::accept_rend() const{
  return accept_iterator(this, -1, false);
}


AliParticleContainer::accept_iterator::accept_iterator(
    const AliParticleContainer *cont,
    int start,
    bool forward):
    fkContainer(cont),
    fAcceptIndices(),
    fCurrentPos(start),
    fForward(forward)
{
  fAcceptIndices.Set(fkContainer->GetNAcceptedParticles());
  UInt_t rejection = 0;
  Int_t mappos = 0;
  for(int i = 0; i < fkContainer->GetNEntries(); i++){
    if(fkContainer->AcceptParticle(i, rejection)){
      fAcceptIndices[mappos++] = i;
    }
  }
}

AliParticleContainer::accept_iterator::accept_iterator(const accept_iterator &other):
    fkContainer(other.fkContainer),
    fAcceptIndices(other.fAcceptIndices),
    fCurrentPos(other.fCurrentPos),
    fForward(other.fForward)
{
}

AliParticleContainer::accept_iterator &AliParticleContainer::accept_iterator::operator=(const AliParticleContainer::accept_iterator &other){
  if(this != &other){
    fkContainer = other.fkContainer;
    fAcceptIndices = other.fAcceptIndices;
    fCurrentPos = other.fCurrentPos;
    fForward = other.fForward;
  }
  return *this;
}

bool AliParticleContainer::accept_iterator::operator!=(const AliParticleContainer::accept_iterator &other) const{
  return fCurrentPos != other.fCurrentPos;
}

AliParticleContainer::accept_iterator &AliParticleContainer::accept_iterator::operator++(){
  if(fForward){
    fCurrentPos++;
  } else {
    fCurrentPos--;
  }
  return *this;
}

AliParticleContainer::accept_iterator AliParticleContainer::accept_iterator::operator++(int){
  AliParticleContainer::accept_iterator result(*this);
  this->operator++();
  return result;
}

AliParticleContainer::accept_iterator &AliParticleContainer::accept_iterator::operator--(){
  if(fForward){
    fCurrentPos--;
  } else {
    fCurrentPos++;
  }
  return *this;
}

AliParticleContainer::accept_iterator AliParticleContainer::accept_iterator::operator--(int){
  AliParticleContainer::accept_iterator result(*this);
  this->operator++();
  return result;
}

AliVParticle *AliParticleContainer::accept_iterator::operator*() const{
  if(fCurrentPos < 0 || fCurrentPos >= fAcceptIndices.GetSize()){
    return NULL;
  }
  return fkContainer->GetParticle(fAcceptIndices[fCurrentPos]);
}


AliParticleContainer::all_iterator::all_iterator(
    const AliParticleContainer *cont,
    int startpos,
    bool forward):
    fkContainer(cont),
    fCurrentPos(startpos),
    fForward(forward)
{
}

AliParticleContainer::all_iterator::all_iterator(const all_iterator &other):
    fkContainer(other.fkContainer),
    fCurrentPos(other.fCurrentPos),
    fForward(other.fForward)
{
}

AliParticleContainer::all_iterator &AliParticleContainer::all_iterator::operator=(const AliParticleContainer::all_iterator &other){
  if(&other != this){
    fkContainer = other.fkContainer;
    fCurrentPos = other.fCurrentPos;
    fForward = other.fForward;
  }
  return *this;
}

bool AliParticleContainer::all_iterator::operator!=(const AliParticleContainer::all_iterator &other) const{
  return other.fCurrentPos != fCurrentPos;
}

AliParticleContainer::all_iterator &AliParticleContainer::all_iterator::operator++(){
  if(fForward){
    fCurrentPos++;
  } else {
    fCurrentPos--;
  }
  return *this;
}

AliParticleContainer::all_iterator AliParticleContainer::all_iterator::operator++(int){
  all_iterator tmp(*this);
  operator++();
  return tmp;
}

AliParticleContainer::all_iterator &AliParticleContainer::all_iterator::operator--(){
  if(fForward){
    fCurrentPos--;
  } else {
    fCurrentPos++;
  }
  return *this;
}

AliParticleContainer::all_iterator AliParticleContainer::all_iterator::operator--(int){
  all_iterator tmp(*this);
  operator--();
  return tmp;
}

AliVParticle *AliParticleContainer::all_iterator::operator*() const{
  if(fCurrentPos >= 0 && fCurrentPos < fkContainer->GetNParticles())
    return fkContainer->GetParticle(fCurrentPos);
  return NULL;
}

/******************************************
 * Unit tests                             *
 ******************************************/

int TestParticleContainerIterator(const AliParticleContainer *const cont, int iteratorType, bool verbose){
  std::vector<AliVParticle *> reference, variation;
  AliVParticle *test = NULL;
  for(int iclust = 0; iclust < cont->GetNParticles(); iclust++){
    test = cont->GetParticle(iclust);
    if(!iteratorType){
      UInt_t rejectionReason = 0;
      if(!cont->AcceptParticle(test, rejectionReason)) continue;
    }
    reference.push_back(test);
  }

  if(!iteratorType){
    // test accept iterator
    for(AliParticleContainer::accept_iterator iter = cont->accept_begin(); iter != cont->accept_end(); ++iter){
      variation.push_back(*iter);
    }
  } else {
    // test all iterator
    for(AliParticleContainer::all_iterator iter = cont->begin(); iter != cont->end(); ++iter){
      variation.push_back(*iter);
    }
  }

  if(verbose) {
    // Printing cluster addresses:
    if(reference.size() < 30){
      std::cout << "Paritcles in reference container: " << std::endl;
      std::cout << "===========================================" << std::endl;
      for(std::vector<AliVParticle *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
        std::cout << "Address: " << *refit << std::endl;
      }
    }
    if(variation.size() < 30){
      std::cout << "Paritcles in test container: " << std::endl;
      std::cout << "===========================================" << std::endl;
      for(std::vector<AliVParticle *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
        std::cout << "Address: " << *varit << std::endl;
      }
    }
  }

  int testresult = 0;
  // compare distributions: all particles in one vector needs to be found in the other vector and vice versa
  bool failure = false;
  // first test: cleck if all particles are found by the iterator
  for(std::vector<AliVParticle *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
    if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
      if(verbose)
        std::cout << "Could not find particle with address " << *clit << " in test container" << std::endl;
      failure = true;
      break;
    }
  }
  if(!failure){
    // second test: check if there are no extra particles found
    for(std::vector<AliVParticle *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
      if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
        if(verbose)
          std::cout << "Could not find particle with address " << *clit << " in reference container" << std::endl;
        failure = true;
        break;
      }
    }
    if(failure) testresult = 2;
  } else testresult = 1;
  if(verbose) {
    std::cout << "Unit test particle container, iterator type " << iteratorType << std::endl;
    std::cout << "Number of expected particles: " << reference.size() << ", number of found particles: " << variation.size() << std::endl;
    std::cout << "Test result: " << testresult << std::endl;
    std::cout << "-----------------------------------------------------------------------" << std::endl;
  }
  return testresult;
}