[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskEmcalJetV0CF.cxx

#include <THnSparse.h>
#include <TClonesArray.h>
#include <TParticle.h>

#include "AliAODEvent.h"
#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliMCParticle.h"
#include "AliAODMCParticle.h"

#include "AliEmcalJet.h"
#include "AliJetContainer.h"
#include "AliParticleContainer.h"
#include "AliClusterContainer.h"

#include "AliPicoV0MC.h"
#include "AliAnalysisTaskEmcalJetV0CF.h"

ClassImp(AliAnalysisTaskEmcalJetV0CF)

//_____________________________________________________________________________
AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF() :
AliAnalysisTaskEmcalJet(),
fKaCutNS(0.),
fLaCutNS(0.),
fV0CutMinEta(0.),
fV0CutMaxEta(0.),
fEventAOD(0),
fEventESD(0),
fCentInfo(0),
fJetsContRD(0),
fTracksContRD(0),
fCaloClustersContRD(0),
fJetsContMC(0),
fTracksContMC(0),
fV0s(0),
fHistoKshortInvM(0),
fHistoLambdaInvM(0),
fHistoAntiLaInvM(0),
fListUserOutputs(0)
{
//
//  AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF
//
}

//_____________________________________________________________________________
AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF(const char *name, Bool_t bHistos) :
AliAnalysisTaskEmcalJet(name,bHistos),
fKaCutNS(6.),
fLaCutNS(6.),
fV0CutMinEta(-10.),
fV0CutMaxEta(10.),
fEventAOD(0),
fEventESD(0),
fCentInfo(0),
fJetsContRD(0),
fTracksContRD(0),
fCaloClustersContRD(0),
fJetsContMC(0),
fTracksContMC(0),
fV0s(0),
fHistoKshortInvM(0),
fHistoLambdaInvM(0),
fHistoAntiLaInvM(0),
fListUserOutputs(0)
{
//
//  AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF
//

  AliAnalysisTaskEmcal::fGeneralHistograms = bHistos;

  DefineOutput(2, TList::Class());
}

//_____________________________________________________________________________
AliAnalysisTaskEmcalJetV0CF::~AliAnalysisTaskEmcalJetV0CF()
{
//
//  AliAnalysisTaskEmcalJetV0CF::~AliAnalysisTaskEmcalJetV0CF
//

  if (fEventAOD) { delete fEventAOD; fEventAOD = 0; }
  if (fEventESD) { delete fEventESD; fEventESD = 0; }
  if (fCentInfo) { delete fCentInfo; fCentInfo = 0; }

  if (fJetsContRD)         { delete fJetsContRD;         fJetsContRD         = 0; }
  if (fTracksContRD)       { delete fTracksContRD;       fTracksContRD       = 0; }
  if (fCaloClustersContRD) { delete fCaloClustersContRD; fCaloClustersContRD = 0; }

  if (fJetsContMC)   { delete fJetsContMC;   fJetsContMC   = 0; }
  if (fTracksContMC) { delete fTracksContMC; fTracksContMC = 0; }

  if (fV0s) { delete fV0s; fV0s = 0; }

  if (fHistoKshortInvM) { delete fHistoKshortInvM; fHistoKshortInvM = 0; }
  if (fHistoLambdaInvM) { delete fHistoLambdaInvM; fHistoLambdaInvM = 0; }
  if (fHistoAntiLaInvM) { delete fHistoAntiLaInvM; fHistoAntiLaInvM = 0; }

  if (fListUserOutputs) { delete fListUserOutputs; fListUserOutputs = 0; }
}

//_____________________________________________________________________________
void AliAnalysisTaskEmcalJetV0CF::Init()
{
//
//  AliAnalysisTaskEmcalJetV0CF::Init
//

  return;
}

//_____________________________________________________________________________
void AliAnalysisTaskEmcalJetV0CF::UserCreateOutputObjects()
{
//
//  AliAnalysisTaskEmcalJetV0CF::UserCreateOutputObjects
//

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
//=============================================================================

  fJetsContRD = GetJetContainer(0);
  if (fJetsContRD) {
    fTracksContRD       = fJetsContRD->GetParticleContainer();
    fCaloClustersContRD = fJetsContRD->GetClusterContainer();
  }

  fJetsContMC = GetJetContainer(1);
  if (fJetsContMC) fTracksContMC = fJetsContMC->GetParticleContainer();
//=============================================================================

  fListUserOutputs = new TList();
  fListUserOutputs->SetOwner();
  CreateUserOutputHistograms();
  PostData(2, fListUserOutputs);
  return;
}

//_____________________________________________________________________________
void AliAnalysisTaskEmcalJetV0CF::Terminate(Option_t *opt)
{
//
//  AliAnalysisTaskEmcalJetV0CF::Terminate
//

  AliAnalysisTaskEmcalJet::Terminate(opt);

  return;
}

//_____________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetV0CF::Run()
{
//
//  AliAnalysisTaskEmcalJetV0CF::Run
//

  if (!AliAnalysisTaskEmcalJet::Run()) return kFALSE;

  return kTRUE;
}

//_____________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetV0CF::RetrieveEventObjects()
{
//
//  AliAnalysisTaskEmcalJetV0CF::RetrieveEventObjects
//

  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) return kFALSE;

  fCentInfo = InputEvent()->GetCentrality(); if (!fCentInfo) return kFALSE;

  fEventAOD = dynamic_cast<AliAODEvent*>(InputEvent());
  fEventESD = dynamic_cast<AliESDEvent*>(InputEvent());
  if ((!fEventAOD) && (!fEventESD)) return kFALSE;

  return kTRUE;
}

//_____________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetV0CF::IsEventSelected()
{
//
//  AliAnalysisTaskEmcalJetV0CF::IsEventSelected
//

  if (!AliAnalysisTaskEmcalJet::IsEventSelected()) return kFALSE;

  return kTRUE;
}

//_____________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetV0CF::FillHistograms()
{
//
//  AliAnalysisTaskEmcalJetV0CF::FillHistograms
//

  if (!AliAnalysisTaskEmcalJet::FillHistograms()) return kFALSE;

  if (FillRecoInfo()) return kFALSE;
  if (FillKineInfo()) return kFALSE;

  return kTRUE;
}

//_____________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetV0CF::FillGeneralHistograms()
{
//
//  AliAnalysisTaskEmcalJetV0CF::FillGeneralHistograms
//

  if (!AliAnalysisTaskEmcalJet::FillGeneralHistograms()) return kFALSE;

  return kTRUE;
}

//_____________________________________________________________________________
void AliAnalysisTaskEmcalJetV0CF::ExecOnce()
{
//
//  AliAnalysisTaskEmcalJetV0CF::ExecOnce
//

  AliAnalysisTaskEmcalJet::ExecOnce();

  if (!fInitialized) return;

  if (!fV0s) {
    fV0s = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoV0s"));

    if (!fV0s) {
      AliError(Form("%s: Could not retrieve V0 %s!", GetName(), "PicoV0s"));
      fInitialized = kFALSE;
      return;
    }
  }

  return;
}

//_____________________________________________________________________________
void AliAnalysisTaskEmcalJetV0CF::CreateUserOutputHistograms()
{
//
//  AliAnalysisTaskEmcalJetV0CF::CreateUserOutputHistograms
//

  if (!fListUserOutputs) return;

  Bool_t bStatusTmpH = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  const Int_t nV0 = 8; // 0: particle type
                       //    ==0, Kshort
                       //    ==1, Lambda
                       //    ==2, AntiLa
                       // 1: Jet pT bin
                       //    == 0.5, jet pT>10.
                       //    == 1.5, jet pT>15.
                       //    == 2.5, jet pT>20.
                       //    == 3.5, jet pT>25.
                       // 2: V0M
                       // 3: V0A
                       // 4: CL1
                       // 5: ZNA
                       // 6: eta
                       // 7: Pt
  const Int_t    nV0Bin[nV0] = {  3,   4,  210,  210,  210,  210, 100, 1000  };
  const Double_t dV0Min[nV0] = { -0.5, 0., -10., -10., -10., -10., -5.,   0. };
  const Double_t dV0Max[nV0] = {  2.5, 4., 200., 200., 200., 200.,  5., 100. };

  THnSparseD *hsReco = new THnSparseD("hsReco", "", nV0, nV0Bin, dV0Min, dV0Max); fListUserOutputs->Add(hsReco);
  THnSparseD *hsKine = new THnSparseD("hsKine", "", nV0, nV0Bin, dV0Min, dV0Max); fListUserOutputs->Add(hsKine);

  TH1::AddDirectory(bStatusTmpH);
  return;
}

//_____________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetV0CF::FillRecoInfo()
{
//
//  AliAnalysisTaskEmcalJetV0CF::FillRecoInfo
//

  if (!fV0s) return kTRUE;
//=============================================================================

  Double_t dV0M = fCentInfo->GetCentralityPercentile("V0M");
  Double_t dV0A = fCentInfo->GetCentralityPercentile("V0A");
  Double_t dCL1 = fCentInfo->GetCentralityPercentile("CL1");
  Double_t dZNA = fCentInfo->GetCentralityPercentile("ZNA");

  THnSparseD *hs = dynamic_cast<THnSparseD*>(fListUserOutputs->FindObject("hsReco"));
  if (!hs)
    return kTRUE;  // should not happen; make Coverity happen
//=============================================================================

  AliPicoV0MC *pV0 = 0;
  for (Int_t i=0; i<fV0s->GetEntriesFast(); i++) {
    pV0 = static_cast<AliPicoV0MC*>(fV0s->At(i)); if (!pV0) continue;

    if (!pV0->IsV0PhysicalPrimary()) { pV0 = 0; continue; }
    if (!pV0->IsV0InEtaAcc(fV0CutMinEta,fV0CutMaxEta)) { pV0 = 0; continue; }

    TVector3 vV0 = pV0->KineMC().Vect();
    Double_t dPt = vV0.Pt();
    Double_t dVar[8];

    dVar[0] = -1.;
    if (pV0->IsKshort()) {
      Int_t k = fHistoKshortInvM->FindBin(dPt); if (k<=0) { pV0 = 0; continue; }

      Double_t dMean  = fHistoKshortInvM->GetBinContent(k);
      Double_t dSigma = fHistoKshortInvM->GetBinError(k);

      Double_t dUpperL = dMean - (fKaCutNS * dSigma);
      Double_t dLowerR = dMean + (fKaCutNS * dSigma);

      Double_t dInvM = pV0->KineKshort().M();
      if ((dInvM<dUpperL) || (dInvM>=dLowerR)) { pV0 = 0; continue; }
      
      dVar[0] = 0.;
    }

    if (pV0->IsLambda()) {
      Int_t k = fHistoLambdaInvM->FindBin(dPt); if (k<=0) { pV0 = 0; continue; }

      Double_t dMean  = fHistoLambdaInvM->GetBinContent(k);
      Double_t dSigma = fHistoLambdaInvM->GetBinError(k);

      Double_t dUpperL = dMean - (fLaCutNS * dSigma);
      Double_t dLowerR = dMean + (fLaCutNS * dSigma);

      Double_t dInvM = pV0->KineLambda().M();
      if ((dInvM<dUpperL) || (dInvM>=dLowerR)) { pV0 = 0; continue; }
      
      dVar[0] = 1.;
    }

    if (pV0->IsAntiLa()) {
      Int_t k = fHistoAntiLaInvM->FindBin(dPt); if (k<=0) { pV0 = 0; continue; }

      Double_t dMean  = fHistoAntiLaInvM->GetBinContent(k);
      Double_t dSigma = fHistoAntiLaInvM->GetBinError(k);

      Double_t dUpperL = dMean - (fLaCutNS * dSigma);
      Double_t dLowerR = dMean + (fLaCutNS * dSigma);

      Double_t dInvM = pV0->KineAntiLa().M();
      if ((dInvM<dUpperL) || (dInvM>=dLowerR)) { pV0 = 0; continue; }

      dVar[0] = 2.;
    }

    if (dVar[0]<-0.5) { pV0 = 0; continue; }

    dVar[1] = -1.;
    if (IsV0InJet(vV0,10.)) dVar[1] = 0.5;
    if (IsV0InJet(vV0,15.)) dVar[1] = 1.5;
    if (IsV0InJet(vV0,20.)) dVar[1] = 2.5;
    if (IsV0InJet(vV0,25.)) dVar[1] = 3.5;
    if (dVar[1]<0.) { pV0 = 0; continue; }

    dVar[2] = dV0M;
    dVar[3] = dV0A;
    dVar[4] = dCL1;
    dVar[5] = dZNA;

    dVar[6] = vV0.Eta();
    dVar[7] = dPt;

    hs->Fill(dVar);

    pV0 = 0;
  }

  return kFALSE;
}

//_____________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetV0CF::FillKineInfo()
{
//
//  AliAnalysisTaskEmcalJetV0CF::FillKineInfo
//

  Double_t dV0M = fCentInfo->GetCentralityPercentile("V0M");
  Double_t dV0A = fCentInfo->GetCentralityPercentile("V0A");
  Double_t dCL1 = fCentInfo->GetCentralityPercentile("CL1");
  Double_t dZNA = fCentInfo->GetCentralityPercentile("ZNA");

  THnSparseD *hs = dynamic_cast<THnSparseD*>(fListUserOutputs->FindObject("hsKine"));
  if (!hs) {
    return kTRUE; // Should not happen; make Coverity happy
  }
//=============================================================================

  AliStack *pStack = 0;
  if (fEventESD) { pStack = MCEvent()->Stack(); if (!pStack) return kTRUE; }
//=============================================================================

  TParticle        *pESD = 0;
  AliAODMCParticle *pAOD = 0;
  for (Int_t i=0; i<MCEvent()->GetNumberOfTracks(); i++) {
    if (fEventAOD) { pAOD = (AliAODMCParticle*)MCEvent()->GetTrack(i);              if (!pAOD) continue; }
    if (fEventESD) { pESD =   ((AliMCParticle*)MCEvent()->GetTrack(i))->Particle(); if (!pESD) continue; }

    Bool_t bPhy = kFALSE;
    if (pAOD) bPhy =   pAOD->IsPhysicalPrimary();
    if (pESD) bPhy = pStack->IsPhysicalPrimary(i);
    if (!bPhy) { pAOD=0; pESD=0; continue; }

    Double_t  dEta = 0.;
    if (pAOD) dEta = pAOD->Eta();
    if (pESD) dEta = pESD->Eta();
    if ((dEta<fV0CutMinEta) || (dEta>=fV0CutMaxEta)) { pAOD=0; pESD=0; continue; }

    Int_t id = 0;
    if (pAOD) id = pAOD->GetPdgCode();
    if (pESD) id = pESD->GetPdgCode();

    Double_t dVar[8]; dVar[0] = -1.;
    if (id== 310 ) dVar[0] = 0.;
    if (id== 3122) dVar[0] = 1.;
    if (id==-3122) dVar[0] = 2.;
    if (dVar[0]<-0.5) { pAOD=0; pESD=0; continue; }

    TVector3 vV0;
    if (pAOD) vV0.SetXYZ(pAOD->Px(), pAOD->Py(), pAOD->Pz());
    if (pESD) vV0.SetXYZ(pESD->Px(), pESD->Py(), pESD->Pz());

    dVar[1] = -1.;
    if (IsV0InJet(vV0,10.)) dVar[1] = 0.5;
    if (IsV0InJet(vV0,15.)) dVar[1] = 1.5;
    if (IsV0InJet(vV0,20.)) dVar[1] = 2.5;
    if (IsV0InJet(vV0,25.)) dVar[1] = 3.5;
    if (dVar[1]<0.) { pAOD=0; pESD=0; continue; }

    dVar[2] = dV0M;
    dVar[3] = dV0A;
    dVar[4] = dCL1;
    dVar[5] = dZNA;
    dVar[6] = dEta;
    if (pAOD) dVar[7] = pAOD->Pt();
    if (pESD) dVar[7] = pESD->Pt();
    hs->Fill(dVar);

    pAOD = 0;
    pESD = 0;
  }

  return kFALSE;
}

//_____________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetV0CF::IsV0InJet(TVector3 vV0, Double_t dJetPtMin)
{
//
//  AliAnalysisTaskEmcalJetV0CF::IsV0InJet
//

  if (!fJetsContRD) return kFALSE;
//=============================================================================

  TVector3 vJet;
  Double_t dJetRadius = fJetsContRD->GetJetRadius();
  AliEmcalJet *pJet = fJetsContRD->GetNextAcceptJet(0); while (pJet) {
    Double_t dPt = fJetsContRD->GetJetPtCorr(fJetsContRD->GetCurrentID());
    if (dPt<dJetPtMin) { pJet = fJetsContRD->GetNextAcceptJet(); continue; }

    vJet.SetPtEtaPhi(dPt, pJet->Eta(), pJet->Phi());
    if (vJet.DeltaR(vV0)<dJetRadius) return kTRUE;
    pJet = fJetsContRD->GetNextAcceptJet();
  }

  return kFALSE;
}
Commit	Line	Data
b254f323	1	#include <THnSparse.h>
	2	#include <TClonesArray.h>
	3	#include <TParticle.h>
	4
	5	#include "AliAODEvent.h"
	6	#include "AliESDEvent.h"
	7	#include "AliMCEvent.h"
	8	#include "AliStack.h"
	9	#include "AliMCParticle.h"
	10	#include "AliAODMCParticle.h"
	11
	12	#include "AliEmcalJet.h"
	13	#include "AliJetContainer.h"
	14	#include "AliParticleContainer.h"
	15	#include "AliClusterContainer.h"
	16
	17	#include "AliPicoV0MC.h"
	18	#include "AliAnalysisTaskEmcalJetV0CF.h"
	19
	20	ClassImp(AliAnalysisTaskEmcalJetV0CF)
	21
	22	//_____________________________________________________________________________
	23	AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF() :
	24	AliAnalysisTaskEmcalJet(),
	25	fKaCutNS(0.),
	26	fLaCutNS(0.),
	27	fV0CutMinEta(0.),
	28	fV0CutMaxEta(0.),
	29	fEventAOD(0),
	30	fEventESD(0),
	31	fCentInfo(0),
	32	fJetsContRD(0),
	33	fTracksContRD(0),
	34	fCaloClustersContRD(0),
	35	fJetsContMC(0),
	36	fTracksContMC(0),
	37	fV0s(0),
	38	fHistoKshortInvM(0),
	39	fHistoLambdaInvM(0),
	40	fHistoAntiLaInvM(0),
	41	fListUserOutputs(0)
	42	{
	43	//
	44	// AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF
	45	//
	46	}
	47
	48	//_____________________________________________________________________________
	49	AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF(const char *name, Bool_t bHistos) :
	50	AliAnalysisTaskEmcalJet(name,bHistos),
	51	fKaCutNS(6.),
	52	fLaCutNS(6.),
	53	fV0CutMinEta(-10.),
	54	fV0CutMaxEta(10.),
	55	fEventAOD(0),
	56	fEventESD(0),
	57	fCentInfo(0),
	58	fJetsContRD(0),
	59	fTracksContRD(0),
	60	fCaloClustersContRD(0),
	61	fJetsContMC(0),
	62	fTracksContMC(0),
	63	fV0s(0),
	64	fHistoKshortInvM(0),
65	fHistoLambdaInvM(0),
66	fHistoAntiLaInvM(0),
67	fListUserOutputs(0)
68	{
69	//
70	// AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF
71	//
72
73	AliAnalysisTaskEmcal::fGeneralHistograms = bHistos;
74
75	DefineOutput(2, TList::Class());
76	}
77
78	//_____________________________________________________________________________
79	AliAnalysisTaskEmcalJetV0CF::~AliAnalysisTaskEmcalJetV0CF()
80	{
81	//
82	// AliAnalysisTaskEmcalJetV0CF::~AliAnalysisTaskEmcalJetV0CF
83	//
84
85	if (fEventAOD) { delete fEventAOD; fEventAOD = 0; }
86	if (fEventESD) { delete fEventESD; fEventESD = 0; }
87	if (fCentInfo) { delete fCentInfo; fCentInfo = 0; }
88
89	if (fJetsContRD) { delete fJetsContRD; fJetsContRD = 0; }
90	if (fTracksContRD) { delete fTracksContRD; fTracksContRD = 0; }
91	if (fCaloClustersContRD) { delete fCaloClustersContRD; fCaloClustersContRD = 0; }
92
93	if (fJetsContMC) { delete fJetsContMC; fJetsContMC = 0; }
94	if (fTracksContMC) { delete fTracksContMC; fTracksContMC = 0; }
95
96	if (fV0s) { delete fV0s; fV0s = 0; }
97
98	if (fHistoKshortInvM) { delete fHistoKshortInvM; fHistoKshortInvM = 0; }
99	if (fHistoLambdaInvM) { delete fHistoLambdaInvM; fHistoLambdaInvM = 0; }
100	if (fHistoAntiLaInvM) { delete fHistoAntiLaInvM; fHistoAntiLaInvM = 0; }
101
102	if (fListUserOutputs) { delete fListUserOutputs; fListUserOutputs = 0; }
103	}
104
105	//_____________________________________________________________________________
106	void AliAnalysisTaskEmcalJetV0CF::Init()
107	{
108	//
109	// AliAnalysisTaskEmcalJetV0CF::Init
110	//
111
112	return;
113	}
114
115	//_____________________________________________________________________________
116	void AliAnalysisTaskEmcalJetV0CF::UserCreateOutputObjects()
117	{
118	//
119	// AliAnalysisTaskEmcalJetV0CF::UserCreateOutputObjects
120	//
121
122	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
123	//=============================================================================
124
125	fJetsContRD = GetJetContainer(0);
126	if (fJetsContRD) {
127	fTracksContRD = fJetsContRD->GetParticleContainer();
128	fCaloClustersContRD = fJetsContRD->GetClusterContainer();
129	}
130
131	fJetsContMC = GetJetContainer(1);
132	if (fJetsContMC) fTracksContMC = fJetsContMC->GetParticleContainer();
133	//=============================================================================
134
135	fListUserOutputs = new TList();
136	fListUserOutputs->SetOwner();
137	CreateUserOutputHistograms();
138	PostData(2, fListUserOutputs);
139	return;
140	}
141
142	//_____________________________________________________________________________
143	void AliAnalysisTaskEmcalJetV0CF::Terminate(Option_t *opt)
144	{
145	//
146	// AliAnalysisTaskEmcalJetV0CF::Terminate
147	//
148
149	AliAnalysisTaskEmcalJet::Terminate(opt);
150
151	return;
152	}
153
154	//_____________________________________________________________________________
155	Bool_t AliAnalysisTaskEmcalJetV0CF::Run()
156	{
157	//
158	// AliAnalysisTaskEmcalJetV0CF::Run
159	//
160
161	if (!AliAnalysisTaskEmcalJet::Run()) return kFALSE;
162
163	return kTRUE;
164	}
165
166	//_____________________________________________________________________________
167	Bool_t AliAnalysisTaskEmcalJetV0CF::RetrieveEventObjects()
168	{
169	//
170	// AliAnalysisTaskEmcalJetV0CF::RetrieveEventObjects
171	//
172
173	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) return kFALSE;
174
175	fCentInfo = InputEvent()->GetCentrality(); if (!fCentInfo) return kFALSE;
176
177	fEventAOD = dynamic_cast<AliAODEvent*>(InputEvent());
178	fEventESD = dynamic_cast<AliESDEvent*>(InputEvent());
179	if ((!fEventAOD) && (!fEventESD)) return kFALSE;
180
181	return kTRUE;
182	}
183
184	//_____________________________________________________________________________
185	Bool_t AliAnalysisTaskEmcalJetV0CF::IsEventSelected()
186	{
187	//
188	// AliAnalysisTaskEmcalJetV0CF::IsEventSelected
189	//
190
191	if (!AliAnalysisTaskEmcalJet::IsEventSelected()) return kFALSE;
192
193	return kTRUE;
194	}
195
196	//_____________________________________________________________________________
197	Bool_t AliAnalysisTaskEmcalJetV0CF::FillHistograms()
198	{
199	//
200	// AliAnalysisTaskEmcalJetV0CF::FillHistograms
201	//
202
203	if (!AliAnalysisTaskEmcalJet::FillHistograms()) return kFALSE;
204
205	if (FillRecoInfo()) return kFALSE;
206	if (FillKineInfo()) return kFALSE;
207
208	return kTRUE;
209	}
210
211	//_____________________________________________________________________________
212	Bool_t AliAnalysisTaskEmcalJetV0CF::FillGeneralHistograms()
213	{
214	//
215	// AliAnalysisTaskEmcalJetV0CF::FillGeneralHistograms
216	//
217
218	if (!AliAnalysisTaskEmcalJet::FillGeneralHistograms()) return kFALSE;
219
220	return kTRUE;
221	}
222
223	//_____________________________________________________________________________
224	void AliAnalysisTaskEmcalJetV0CF::ExecOnce()
225	{
226	//
227	// AliAnalysisTaskEmcalJetV0CF::ExecOnce
228	//
229
230	AliAnalysisTaskEmcalJet::ExecOnce();
231
232	if (!fInitialized) return;
233
234	if (!fV0s) {
235	fV0s = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoV0s"));
236
237	if (!fV0s) {
238	AliError(Form("%s: Could not retrieve V0 %s!", GetName(), "PicoV0s"));
239	fInitialized = kFALSE;
240	return;
241	}
242	}
243
244	return;
245	}
246
247	//_____________________________________________________________________________
248	void AliAnalysisTaskEmcalJetV0CF::CreateUserOutputHistograms()
249	{
250	//
251	// AliAnalysisTaskEmcalJetV0CF::CreateUserOutputHistograms
252	//
253
254	if (!fListUserOutputs) return;
255
256	Bool_t bStatusTmpH = TH1::AddDirectoryStatus();
257	TH1::AddDirectory(kFALSE);
258
259	const Int_t nV0 = 8; // 0: particle type
260	// ==0, Kshort
261	// ==1, Lambda
262	// ==2, AntiLa
263	// 1: Jet pT bin
264	// == 0.5, jet pT>10.
265	// == 1.5, jet pT>15.
266	// == 2.5, jet pT>20.
267	// == 3.5, jet pT>25.
268	// 2: V0M
269	// 3: V0A
270	// 4: CL1
271	// 5: ZNA
272	// 6: eta
273	// 7: Pt
274	const Int_t nV0Bin[nV0] = { 3, 4, 210, 210, 210, 210, 100, 1000 };
275	const Double_t dV0Min[nV0] = { -0.5, 0., -10., -10., -10., -10., -5., 0. };
276	const Double_t dV0Max[nV0] = { 2.5, 4., 200., 200., 200., 200., 5., 100. };
277
278	THnSparseD *hsReco = new THnSparseD("hsReco", "", nV0, nV0Bin, dV0Min, dV0Max); fListUserOutputs->Add(hsReco);
279	THnSparseD *hsKine = new THnSparseD("hsKine", "", nV0, nV0Bin, dV0Min, dV0Max); fListUserOutputs->Add(hsKine);
280
281	TH1::AddDirectory(bStatusTmpH);
282	return;
283	}
284
285	//_____________________________________________________________________________
286	Bool_t AliAnalysisTaskEmcalJetV0CF::FillRecoInfo()
287	{
288	//
289	// AliAnalysisTaskEmcalJetV0CF::FillRecoInfo
290	//
291
292	if (!fV0s) return kTRUE;
293	//=============================================================================
294
295	Double_t dV0M = fCentInfo->GetCentralityPercentile("V0M");
296	Double_t dV0A = fCentInfo->GetCentralityPercentile("V0A");
297	Double_t dCL1 = fCentInfo->GetCentralityPercentile("CL1");
298	Double_t dZNA = fCentInfo->GetCentralityPercentile("ZNA");
299
300	THnSparseD hs = dynamic_cast<THnSparseD>(fListUserOutputs->FindObject("hsReco"));
81fe377a	301	if (!hs)
81fe377a	302	return kTRUE; // should not happen; make Coverity happen
b254f323	303	//=============================================================================
	304
	305	AliPicoV0MC *pV0 = 0;
	306	for (Int_t i=0; i<fV0s->GetEntriesFast(); i++) {
	307	pV0 = static_cast<AliPicoV0MC*>(fV0s->At(i)); if (!pV0) continue;
	308
	309	if (!pV0->IsV0PhysicalPrimary()) { pV0 = 0; continue; }
	310	if (!pV0->IsV0InEtaAcc(fV0CutMinEta,fV0CutMaxEta)) { pV0 = 0; continue; }
	311
	312	TVector3 vV0 = pV0->KineMC().Vect();
	313	Double_t dPt = vV0.Pt();
	314	Double_t dVar[8];
	315
	316	dVar[0] = -1.;
	317	if (pV0->IsKshort()) {
	318	Int_t k = fHistoKshortInvM->FindBin(dPt); if (k<=0) { pV0 = 0; continue; }
	319
	320	Double_t dMean = fHistoKshortInvM->GetBinContent(k);
	321	Double_t dSigma = fHistoKshortInvM->GetBinError(k);
	322
	323	Double_t dUpperL = dMean - (fKaCutNS * dSigma);
	324	Double_t dLowerR = dMean + (fKaCutNS * dSigma);
	325
	326	Double_t dInvM = pV0->KineKshort().M();
	327	if ((dInvM<dUpperL) \|\| (dInvM>=dLowerR)) { pV0 = 0; continue; }
	328
	329	dVar[0] = 0.;
	330	}
	331
	332	if (pV0->IsLambda()) {
	333	Int_t k = fHistoLambdaInvM->FindBin(dPt); if (k<=0) { pV0 = 0; continue; }
	334
	335	Double_t dMean = fHistoLambdaInvM->GetBinContent(k);
	336	Double_t dSigma = fHistoLambdaInvM->GetBinError(k);
	337
	338	Double_t dUpperL = dMean - (fLaCutNS * dSigma);
	339	Double_t dLowerR = dMean + (fLaCutNS * dSigma);
	340
	341	Double_t dInvM = pV0->KineLambda().M();
	342	if ((dInvM<dUpperL) \|\| (dInvM>=dLowerR)) { pV0 = 0; continue; }
	343
	344	dVar[0] = 1.;
	345	}
	346
	347	if (pV0->IsAntiLa()) {
	348	Int_t k = fHistoAntiLaInvM->FindBin(dPt); if (k<=0) { pV0 = 0; continue; }
	349
	350	Double_t dMean = fHistoAntiLaInvM->GetBinContent(k);
	351	Double_t dSigma = fHistoAntiLaInvM->GetBinError(k);
	352
	353	Double_t dUpperL = dMean - (fLaCutNS * dSigma);
	354	Double_t dLowerR = dMean + (fLaCutNS * dSigma);
	355
	356	Double_t dInvM = pV0->KineAntiLa().M();
	357	if ((dInvM<dUpperL) \|\| (dInvM>=dLowerR)) { pV0 = 0; continue; }
	358
	359	dVar[0] = 2.;
	360	}
	361
	362	if (dVar[0]<-0.5) { pV0 = 0; continue; }
	363
	364	dVar[1] = -1.;
	365	if (IsV0InJet(vV0,10.)) dVar[1] = 0.5;
	366	if (IsV0InJet(vV0,15.)) dVar[1] = 1.5;
367	if (IsV0InJet(vV0,20.)) dVar[1] = 2.5;
368	if (IsV0InJet(vV0,25.)) dVar[1] = 3.5;
369	if (dVar[1]<0.) { pV0 = 0; continue; }
370
371	dVar[2] = dV0M;
372	dVar[3] = dV0A;
373	dVar[4] = dCL1;
374	dVar[5] = dZNA;
375
376	dVar[6] = vV0.Eta();
377	dVar[7] = dPt;
378
379	hs->Fill(dVar);
380
381	pV0 = 0;
382	}
383
384	return kFALSE;
385	}
386
387	//_____________________________________________________________________________
388	Bool_t AliAnalysisTaskEmcalJetV0CF::FillKineInfo()
389	{
390	//
391	// AliAnalysisTaskEmcalJetV0CF::FillKineInfo
392	//
393
394	Double_t dV0M = fCentInfo->GetCentralityPercentile("V0M");
395	Double_t dV0A = fCentInfo->GetCentralityPercentile("V0A");
396	Double_t dCL1 = fCentInfo->GetCentralityPercentile("CL1");
397	Double_t dZNA = fCentInfo->GetCentralityPercentile("ZNA");
398
399	THnSparseD hs = dynamic_cast<THnSparseD>(fListUserOutputs->FindObject("hsKine"));
81fe377a	400	if (!hs) {
	401	return kTRUE; // Should not happen; make Coverity happy
	402	}
b254f323	403	//=============================================================================
	404
	405	AliStack *pStack = 0;
	406	if (fEventESD) { pStack = MCEvent()->Stack(); if (!pStack) return kTRUE; }
	407	//=============================================================================
	408
	409	TParticle *pESD = 0;
	410	AliAODMCParticle *pAOD = 0;
	411	for (Int_t i=0; i<MCEvent()->GetNumberOfTracks(); i++) {
	412	if (fEventAOD) { pAOD = (AliAODMCParticle*)MCEvent()->GetTrack(i); if (!pAOD) continue; }
	413	if (fEventESD) { pESD = ((AliMCParticle*)MCEvent()->GetTrack(i))->Particle(); if (!pESD) continue; }
	414
	415	Bool_t bPhy = kFALSE;
	416	if (pAOD) bPhy = pAOD->IsPhysicalPrimary();
	417	if (pESD) bPhy = pStack->IsPhysicalPrimary(i);
	418	if (!bPhy) { pAOD=0; pESD=0; continue; }
	419
	420	Double_t dEta = 0.;
	421	if (pAOD) dEta = pAOD->Eta();
	422	if (pESD) dEta = pESD->Eta();
	423	if ((dEta<fV0CutMinEta) \|\| (dEta>=fV0CutMaxEta)) { pAOD=0; pESD=0; continue; }
	424
	425	Int_t id = 0;
	426	if (pAOD) id = pAOD->GetPdgCode();
	427	if (pESD) id = pESD->GetPdgCode();
	428
	429	Double_t dVar[8]; dVar[0] = -1.;
	430	if (id== 310 ) dVar[0] = 0.;
	431	if (id== 3122) dVar[0] = 1.;
	432	if (id==-3122) dVar[0] = 2.;
	433	if (dVar[0]<-0.5) { pAOD=0; pESD=0; continue; }
	434
	435	TVector3 vV0;
	436	if (pAOD) vV0.SetXYZ(pAOD->Px(), pAOD->Py(), pAOD->Pz());
	437	if (pESD) vV0.SetXYZ(pESD->Px(), pESD->Py(), pESD->Pz());
	438
	439	dVar[1] = -1.;
	440	if (IsV0InJet(vV0,10.)) dVar[1] = 0.5;
	441	if (IsV0InJet(vV0,15.)) dVar[1] = 1.5;
	442	if (IsV0InJet(vV0,20.)) dVar[1] = 2.5;
	443	if (IsV0InJet(vV0,25.)) dVar[1] = 3.5;
	444	if (dVar[1]<0.) { pAOD=0; pESD=0; continue; }
	445
	446	dVar[2] = dV0M;
	447	dVar[3] = dV0A;
	448	dVar[4] = dCL1;
	449	dVar[5] = dZNA;
	450	dVar[6] = dEta;
	451	if (pAOD) dVar[7] = pAOD->Pt();
	452	if (pESD) dVar[7] = pESD->Pt();
	453	hs->Fill(dVar);
	454
	455	pAOD = 0;
	456	pESD = 0;
	457	}
	458
	459	return kFALSE;
	460	}
	461
	462	//_____________________________________________________________________________
	463	Bool_t AliAnalysisTaskEmcalJetV0CF::IsV0InJet(TVector3 vV0, Double_t dJetPtMin)
	464	{
	465	//
	466	// AliAnalysisTaskEmcalJetV0CF::IsV0InJet
467	//
468
469	if (!fJetsContRD) return kFALSE;
470	//=============================================================================
471
472	TVector3 vJet;
473	Double_t dJetRadius = fJetsContRD->GetJetRadius();
474	AliEmcalJet *pJet = fJetsContRD->GetNextAcceptJet(0); while (pJet) {
475	Double_t dPt = fJetsContRD->GetJetPtCorr(fJetsContRD->GetCurrentID());
476	if (dPt<dJetPtMin) { pJet = fJetsContRD->GetNextAcceptJet(); continue; }
477
478	vJet.SetPtEtaPhi(dPt, pJet->Eta(), pJet->Phi());
479	if (vJet.DeltaR(vV0)<dJetRadius) return kTRUE;
480	pJet = fJetsContRD->GetNextAcceptJet();
481	}
482
483	return kFALSE;
484	}