AliEMCalTriggerRecTrackAnalysisComponent.cxx

Go to the documentation of this file.

/**************************************************************************
 * Copyright(c) 1998-2014, 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 <map>
#include <string>
#include <vector>

#include <TAxis.h>
#include <TClonesArray.h>
#include <TList.h>
#include <TMath.h>
#include <TString.h>

#include "AliAODMCParticle.h"
#include "AliLog.h"
#include "AliEmcalTriggerPatchInfo.h"
#include "AliMCEvent.h"
#include "AliPicoTrack.h"
#include "AliVCluster.h"
#include "AliVEvent.h"
#include "AliVParticle.h"
#include "AliVTrack.h"

#include "AliEMCalTriggerAnaTriggerDecision.h"
#include "AliEMCalTriggerBinningComponent.h"
#include "AliEMCalTriggerEventData.h"
#include "AliEMCalTriggerKineCuts.h"
#include "AliEMCalPtTaskVTrackSelection.h"
#include "AliEMCalTriggerRecTrackAnalysisComponent.h"
#include "AliEMCalTriggerWeightHandler.h"

ClassImp(EMCalTriggerPtAnalysis::AliEMCalTriggerRecTrackAnalysisComponent)

namespace EMCalTriggerPtAnalysis {

AliEMCalTriggerRecTrackAnalysisComponent::AliEMCalTriggerRecTrackAnalysisComponent() :
  AliEMCalTriggerTracksAnalysisComponent(),
  fTrackSelection(NULL),
  fSwapEta(kFALSE),
  fTriggerMethod(kTriggerString),
  fRequestMCtrue(kFALSE),
  fDoMatchPatches(kFALSE)
{
}

AliEMCalTriggerRecTrackAnalysisComponent::AliEMCalTriggerRecTrackAnalysisComponent(const char *name) :
  AliEMCalTriggerTracksAnalysisComponent(name),
  fTrackSelection(NULL),
  fSwapEta(kFALSE),
  fTriggerMethod(kTriggerString),
  fRequestMCtrue(kFALSE),
  fDoMatchPatches(kFALSE)
{
}

AliEMCalTriggerRecTrackAnalysisComponent::~AliEMCalTriggerRecTrackAnalysisComponent() {
  if(fTrackSelection) delete fTrackSelection;
}

void AliEMCalTriggerRecTrackAnalysisComponent::CreateHistos() {
  AliEMCalTriggerTracksAnalysisComponent::CreateHistos();

  // Create trigger definitions
  std::map<std::string, std::string> triggerCombinations;
  const char *triggernames[11] = {"MinBias", "EMCJHigh", "EMCJLow", "EMCGHigh",
      "EMCGLow", "EMCHighBoth", "EMCHighGammaOnly", "EMCHighJetOnly",
      "EMCLowBoth", "EMCLowGammaOnly", "EMCLowJetOnly"};
  // Define names and titles for different triggers in the histogram container
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[0], "min. bias events"));
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[1], "jet-triggered events (high threshold)"));
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[2], "jet-triggered events (low threshold)"));
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[3], "gamma-triggered events (high threshold)"));
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[4], "gamma-triggered events (low threshold)"));
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[5], "jet and gamma triggered events (high threshold)"));
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[6], "exclusively gamma-triggered events (high threshold)"));
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[7], "exclusively jet-triggered events (high threshold)"));
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[8], "jet and gamma triggered events (low threshold)"));
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[9], "exclusively gamma-triggered events (low threshold)"));
  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[10], "exclusively-triggered events (low threshold)"));

  // Create axis definitions
  const AliEMCalTriggerBinningDimension *ptbinning = fBinning->GetBinning("pt"),
      *etabinning = fBinning->GetBinning("eta"),
      *phibinning = fBinning->GetBinning("phi"),
      *vertexbinning = fBinning->GetBinning("zvertex");

  const TAxis *trackaxes[5] = {
      DefineAxis("pt", ptbinning),
      DefineAxis("eta", etabinning),
      DefineAxis("phi", phibinning),
      DefineAxis("zvertex", vertexbinning),
      DefineAxis("mbtrigger", 2, -0.5, 1.5)
  };

  // Build histograms
  for(std::map<std::string,std::string>::iterator it = triggerCombinations.begin(); it != triggerCombinations.end(); ++it){
    const std::string name = it->first, &title = it->second;
    fHistos->CreateTHnSparse(Form("hTrackHist%s", name.c_str()), Form("Track-based data for %s events", title.c_str()), 5, trackaxes, "s");
    fHistos->CreateTHnSparse(Form("hTrackInAcceptanceHist%s", name.c_str()), Form("Track-based data for %s events for tracks matched to EMCal clusters", title.c_str()), 5, trackaxes, "s");
    fHistos->CreateTHnSparse(Form("hMCTrackHist%s", name.c_str()), Form("Track-based data for %s events with MC kinematics", title.c_str()), 5, trackaxes, "s");
    fHistos->CreateTHnSparse(Form("hMCTrackInAcceptanceHist%s", name.c_str()), Form("Track-based data for %s events for tracks matched to EMCal clusters with MC kinematics", title.c_str()), 5, trackaxes, "s");
    if(fDoMatchPatches){
      fHistos->CreateTHnSparse(Form("hTrackHistPatchMatch%s", name.c_str()), Form("Track-based data for %s events for tracks matched to EMCal patches", title.c_str()), 5, trackaxes, "s");
      fHistos->CreateTHnSparse(Form("hMCTrackHistPatchMatch%s", name.c_str()), Form("Track-based data for %s events for tracks matched to EMCal patches with MC kinematics", title.c_str()), 5, trackaxes, "s");
    }
  }

  // Correlation Matrix
  // 1. gen pt
  // 2. rec pt
  // 3. rec eta
  // 4. rec phi
  const TAxis *corraxes[4] = {
    DefineAxis("ptgen", ptbinning),
    DefineAxis("ptrec", ptbinning),
    DefineAxis("eta", etabinning),
    DefineAxis("phi", phibinning)
  };
  fHistos->CreateTHnSparse("hTrackPtCorrelation", "Correlation matrix for track pt", 4, corraxes);

  for(int iaxis = 0; iaxis < 5; iaxis++) delete trackaxes[iaxis];
}

void AliEMCalTriggerRecTrackAnalysisComponent::Process(const AliEMCalTriggerEventData* const data) {
  AliDebug(1, Form("Number of matched tracks: %d", data->GetMatchedTrackContainer()->GetEntries()));
  if(fRequestMCtrue && !data->GetMCEvent()) return;

  std::vector<std::string> triggernames;
  this->GetMachingTriggerNames(triggernames, fTriggerMethod);

  AliVTrack *track(NULL);
  const AliVParticle *assocMC(NULL);
  if(!data->GetMatchedTrackContainer()){
    AliError("No container for matched tracks");
    return;
  }
  TIter trackIter(data->GetMatchedTrackContainer());

  double weight = 1.;
  if(fWeightHandler && data->GetMCEvent()){
    weight = fWeightHandler->GetEventWeight(data->GetMCEvent());
  }
  while((track = dynamic_cast<AliVTrack *>(trackIter()))){
    // Apply track selection
    assocMC = NULL;
    if(fKineCuts && !fKineCuts->IsSelected(track)) continue;
    if(fTrackSelection && !fTrackSelection->IsTrackAccepted(track)){
      AliDebug(2, "Track not accepted");
      continue;
    }

    if(fRequestMCtrue){
      if(!(assocMC = IsMCTrueTrack(track, data->GetMCEvent()))) continue; // Not a true track
      this->FillCorrelation(assocMC, track);
    }
    // Try to match the cluster
    Bool_t hasCluster = kFALSE;
    AliVCluster *clust(NULL);
    AliVTrack *testtrack = track;
    AliPicoTrack *pictrack = dynamic_cast<AliPicoTrack *>(track);
    if(pictrack) testtrack = pictrack->GetTrack();
    if(testtrack->GetEMCALcluster() >= 0 && (clust = dynamic_cast<AliVCluster *>(data->GetClusterContainer()->At(testtrack->GetEMCALcluster()))))
      hasCluster = kTRUE;

    // Try to match to EMCAL patches (only in case the event is triggered)
    TList patches;
    if(fDoMatchPatches && triggernames.size()) MatchTriggerPatches(track, data->GetTriggerPatchContainer(), patches);

    // Fill histograms
    for(std::vector<std::string>::iterator name = triggernames.begin(); name != triggernames.end(); name++){
      FillHistogram(Form("hTrackHist%s", name->c_str()), track, NULL, data->GetRecEvent(), kFALSE, weight);
      if(hasCluster) FillHistogram(Form("hTrackInAcceptanceHist%s", name->c_str()), track, NULL, data->GetRecEvent(), kFALSE, weight);
      Bool_t hasPatch = HasMatchedPatchOfType(*name, patches);
      if(fDoMatchPatches && hasPatch) FillHistogram(Form("hTrackHistPatchMatch%s", name->c_str()), track, NULL, data->GetRecEvent(), kFALSE, weight);
      if(assocMC){
        FillHistogram(Form("hMCTrackHist%s", name->c_str()), track, assocMC, data->GetRecEvent(), kTRUE, weight);
        if(hasCluster) FillHistogram(Form("hMCTrackInAcceptanceHist%s", name->c_str()), track, assocMC, data->GetRecEvent(), kTRUE, weight);
        if(fDoMatchPatches && hasPatch) FillHistogram(Form("hMCTrackHistPatchMatch%s", name->c_str()), track, NULL, data->GetRecEvent(), kFALSE, weight);
      }
    }
  }
}

const AliVParticle * AliEMCalTriggerRecTrackAnalysisComponent::IsMCTrueTrack(
    const AliVTrack* const trk, const AliMCEvent* evnt) const {
  int label = TMath::Abs(trk->GetLabel());
  const AliVParticle *mcpart = evnt->GetTrack(label);
  if(!mcpart) return NULL;
  const AliAODMCParticle *aodpart = dynamic_cast<const AliAODMCParticle *>(mcpart);
  if(aodpart){
    if(!aodpart->IsPhysicalPrimary()) return NULL;
  } else {
    if(!evnt->IsPhysicalPrimary(label)) return NULL;
  }
  return mcpart;
}

void AliEMCalTriggerRecTrackAnalysisComponent::FillHistogram(
    const TString& histname, const AliVTrack* const trk,
    const AliVParticle* assocMC, const AliVEvent* const recev,
    Bool_t useMCkine, Double_t weight) {
  if(useMCkine && !assocMC) return;
  double data[5];
  data[0] = useMCkine ? TMath::Abs(assocMC->Pt()) : TMath::Abs(trk->Pt());
  data[1] = (fSwapEta ? -1. : 1.) * (useMCkine ? assocMC->Eta() : trk->Eta());
  data[2] = useMCkine ? assocMC->Phi() : trk->Phi();
  data[3] = recev->GetPrimaryVertex()->GetZ();
  data[4] = fTriggerDecision->IsMinBias();
  fHistos->FillTHnSparse(histname.Data(), data, weight);
}

void AliEMCalTriggerRecTrackAnalysisComponent::FillCorrelation(
    const AliVParticle* const genparticle,
    const AliVParticle* const recparticle, double weight) {
  double data[4] = {TMath::Abs(genparticle->Pt()), TMath::Abs(recparticle->Pt()), recparticle->Eta(), recparticle->Phi()};
  fHistos->FillTHnSparse("hTrackPtCorrelation", data, weight);
}

void AliEMCalTriggerRecTrackAnalysisComponent::MatchTriggerPatches(
    const AliVTrack* rectrack, const TClonesArray* inputpatches,
    TList& outputpatches) const {
  outputpatches.Clear();
  Double_t etaEMCAL = rectrack->GetTrackEtaOnEMCal(),
           phiEMCAL = rectrack->GetTrackPhiOnEMCal();
  for(TIter patchiter = TIter(inputpatches).Begin(); patchiter != TIter::End(); ++patchiter){
    AliEmcalTriggerPatchInfo *recpatch = static_cast<AliEmcalTriggerPatchInfo *>(*patchiter);
    double etamin = TMath::Min(recpatch->GetEtaMin(), recpatch->GetEtaMax()),
        etamax = TMath::Max(recpatch->GetEtaMin(), recpatch->GetEtaMax()),
        phimin = TMath::Min(recpatch->GetPhiMin(), recpatch->GetPhiMax()),
        phimax = TMath::Max(recpatch->GetPhiMin(), recpatch->GetPhiMax());
    if(etaEMCAL >= etamin && etaEMCAL <= etamax && phiEMCAL >= phimin && phiEMCAL <= phimax){
      // Patch accepted
      outputpatches.Add(recpatch);
    }
  }
}

Bool_t AliEMCalTriggerRecTrackAnalysisComponent::HasMatchedPatchOfType(
    TString triggertype, const TList& patches) const {
  Bool_t hasGammaLow = kFALSE,
      hasGammaHigh = kFALSE,
      hasJetLow = kFALSE,
      hasJetHigh = kFALSE;
  for(TIter patchIter = TIter(&patches).Begin(); patchIter != TIter::End(); ++patchIter){
    AliEmcalTriggerPatchInfo *matchedpad = static_cast<AliEmcalTriggerPatchInfo *>(*patchIter);
    if(fTriggerDecision->GetConfiguration()->IsUsingOfflinePatches()){
      if(!matchedpad->IsOfflineSimple()) continue;
      if(matchedpad->IsJetHighSimple()) hasJetHigh = kTRUE;
      if(matchedpad->IsJetLowSimple()) hasJetLow = kTRUE;
      if(matchedpad->IsGammaHighSimple()) hasGammaHigh = kTRUE;
      if(matchedpad->IsGammaLowSimple()) hasGammaLow = kTRUE;
    } else {
      if(matchedpad->IsOfflineSimple()) continue;
      if(matchedpad->IsJetHigh()) hasJetHigh = kTRUE;
      if(matchedpad->IsJetLow()) hasJetLow = kTRUE;
      if(matchedpad->IsGammaHigh()) hasGammaHigh = kTRUE;
      if(matchedpad->IsGammaLow()) hasGammaLow = kTRUE;
    }
  }
  // Evaluate status
  if((triggertype.Contains("EMCJHigh") || triggertype.Contains("EMCHighJetOnly")) && hasJetHigh) return kTRUE;
  if((triggertype.Contains("EMCJLow") || triggertype.Contains("EMCLowJetOnly")) && hasJetLow) return kTRUE;
  if((triggertype.Contains("EMCGHigh") || triggertype.Contains("EMCHighGammaOnly")) && hasGammaHigh) return kTRUE;
  if((triggertype.Contains("EMCGLow") || triggertype.Contains("EMCLowGammaOnly")) && hasGammaLow) return kTRUE;
  if(triggertype.Contains("EMCHighBoth") && (hasJetHigh || hasGammaHigh)) return kTRUE;
  if(triggertype.Contains("EMCLowBoth") && (hasJetLow || hasGammaLow)) return kTRUE;
  return kFALSE;
}

} /* namespace EMCalTriggerPtAnalysis */