[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskdPhi.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *																          *
 * Authors: Svein Lindal                                                  *
 * Version 1.0			                                                  *
 *																          *
 *																          *
 * 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.				  *
 **************************************************************************/

////////////////////////////////////////////////
//--------------------------------------------- 
// Class doing conversion gamma dPhi correlations
// Gamma Conversion analysis
//---------------------------------------------
////////////////////////////////////////////////

#include "AliAnalysisTaskdPhi.h"

#include <TH2I.h>
#include <TList.h>
#include <TChain.h>

#include <AliAnalysisManager.h>
#include <AliInputEventHandler.h>
#include <AliESDInputHandler.h>
#include <AliAODInputHandler.h>
#include <AliAnalysisFilter.h>

#include "AliConversionCuts.h"
#include "AliAODConversionPhoton.h"
#include "AliAODConversionMother.h"
#include "AliAnaConvCorrPhoton.h"
#include "AliAnaConvCorrPion.h"
#include "AliAnaConvIsolation.h"
// Author Svein Lindal <slindal@fys.uio.no>
using namespace std;

ClassImp(AliAnalysisTaskdPhi)


//________________________________________________________________________
AliAnalysisTaskdPhi::AliAnalysisTaskdPhi(const char *name) : AliAnalysisTaskSE(name),
  fHistograms(NULL),
  fHistoGamma(NULL),
  fHistoPion(NULL),
  fDielV0TrackFilter(NULL), 
  fDielV0Filter(NULL),
  fDielPi0Filter(NULL),
  fDielTrackFilter(NULL),
  fV0Filter(NULL),
  fGammas(NULL),
  fPions(NULL),
  hMETracks(NULL), 
  hMEPhotons(NULL), 
  hMEPions(NULL),
  hMEvents(NULL),
  fPhotonCorr(NULL),
  fPionCorr(NULL), 
  fIsoAna(NULL),
  fL1(-1),
  fL2(-1),
  fDeltaAODBranchName("AliAODGammaConversion_gamma"), 
  fAxisPt(),
  fAxisEta(),
  fAxisPhi(),
  fAxisCent(),
  fAxisZ(), 
  fAxisPiM()
{
  //constructor


  fAxisPt.SetNameTitle("PtAxis", "Pt");
  fAxisPt.Set(20, 0, 100);
  fAxisEta.SetNameTitle("EtaAxis", "Eta");
  fAxisEta.Set(160, -0.8, 0.8);
  fAxisPhi.SetNameTitle("PhiAxis", "Phi");
  fAxisPhi.Set(128, 0, TMath::TwoPi());
  fAxisZ.SetNameTitle("ZAxis", "Z");
  fAxisZ.Set(4, -10, 10);
  fAxisCent.SetNameTitle("CentAxis", "Cent");

  Double_t centbins[5] = {0, 10, 30, 60, 100.1};
  fAxisCent.Set(4, centbins);


  Double_t mbins[7] = {0.1, 0.11, 0.12, 0.15, 0.16, 0.18, 0.2};
  fAxisPiM.SetNameTitle("InvMassPi0", "Invariant mass");
  fAxisPiM.Set(6, mbins);

  fGammas = new TObjArray();
  fGammas->SetOwner(kFALSE);

  fPions = new TObjArray();
  fPions->SetOwner(kFALSE);
  
  // Define input and output slots here
  DefineInput(0, TChain::Class());
  DefineInput(1, TClonesArray::Class());
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
  DefineOutput(3, TList::Class());
}


//________________________________________________________________________
AliAnalysisTaskdPhi::~AliAnalysisTaskdPhi(){
  //destructor
  if(fPions)
	delete fPions;
  fPions = NULL;

  if(fGammas)
	delete fGammas;
  fGammas = NULL;
  
  if(fIsoAna)
  	delete fIsoAna;
  fIsoAna = NULL;

  if(fV0Filter)
	delete fV0Filter;
  fV0Filter = NULL;

  if(fHistograms)
	delete fHistograms;
  fHistograms = NULL;

if(fHistoPion)
	delete fHistoPion;
  fHistoPion = NULL;

  if(fHistoGamma)
	delete fHistoGamma;
  fHistoGamma = NULL;

}

///________________________________________________________________________
void AliAnalysisTaskdPhi::SetUpCorrObjects() {
  ///Creat corr obj
  fIsoAna = new AliAnaConvIsolation();


  fPhotonCorr = new TObjArray();
  fPionCorr = new TObjArray();
  
  TList * hPhoton = new TList();
  hPhoton->SetName("hPhotonCorr");
  hPhoton->SetOwner(kTRUE);
  fHistoGamma->Add(hPhoton);

  TList * hPion = new TList();
  hPion->SetName("hPionCorr");
  hPion->SetOwner(kTRUE);
  fHistoPion->Add(hPion);


  for(Int_t iz = 0; iz < fAxisZ.GetNbins(); iz++) {
	TObjArray * photonArray = new TObjArray();
	photonArray->SetOwner(kTRUE);
	fPhotonCorr->AddAt(photonArray, iz);

	TObjArray * pionArray = new TObjArray();
	pionArray->SetOwner(kTRUE);
	fPionCorr->AddAt(pionArray, iz);

	TList * photonList = new TList();
	photonList->SetName(Form("photon_%d", iz));
	photonList->SetOwner(kTRUE);
	hPhoton->AddAt(photonList, iz);

	TList * pionList = new TList();
	pionList->SetName(Form("pion_%d", iz));
	pionList->SetOwner(kTRUE);
	hPion->AddAt(pionList, iz);
	
	for(Int_t ic = 0; ic < fAxisCent.GetNbins(); ic++) {
	  
	  TString nameString = Form("%d_%d", iz, ic);
	  TString titleString = Form("%f < Z < %f ... %f cent %f", 
								 fAxisZ.GetBinLowEdge(iz+1), fAxisZ.GetBinUpEdge(iz+1), 
								 fAxisCent.GetBinLowEdge(ic+1), fAxisCent.GetBinUpEdge(ic+1));


	  AliAnaConvCorrPhoton * photonCorr = new AliAnaConvCorrPhoton(Form("PhotonCorr_%s", nameString.Data()), Form("photon %s", titleString.Data()));
	  photonArray->AddAt(photonCorr, ic);
	  photonCorr->GetAxistPt().Set(fAxisPt.GetNbins(), fAxisPt.GetXbins()->GetArray());
	  photonCorr->GetAxiscPt().Set(fAxisPt.GetNbins(), fAxisPt.GetXbins()->GetArray());
	  photonCorr->CreateHistograms();
	  photonList->Add(photonCorr->GetHistograms());

	  AliAnaConvCorrPion * pionCorr = new AliAnaConvCorrPion(Form("PionCorr_%s", nameString.Data()), Form("pion %s", titleString.Data()));
	  pionArray->AddAt(pionCorr, ic);
	  pionCorr->GetAxistPt().Set(fAxisPt.GetNbins(), fAxisPt.GetXbins()->GetArray());
	  pionCorr->GetAxiscPt().Set(fAxisPt.GetNbins(), fAxisPt.GetXbins()->GetArray());
	  pionCorr->GetAxisM().Set(fAxisPiM.GetNbins(), fAxisPiM.GetXbins()->GetArray());
	  pionCorr->CreateHistograms();
	  pionList->Add(pionCorr->GetHistograms());
	}
  }
}

//________________________________________________________________________
void AliAnalysisTaskdPhi::UserCreateOutputObjects() {
  // Create histograms
  
  fHistograms = new TList();
  fHistograms->SetName("dPhi_histograms");
  fHistograms->SetOwner(kTRUE);

  fHistoGamma = new TList();
  fHistoGamma->SetName("Gamma_histo");
  fHistoGamma->SetOwner(kTRUE);

  fHistoPion = new TList();
  fHistoPion->SetName("Pion_histo");
  fHistoPion->SetOwner(kTRUE);

  
  if(fV0Filter) {
	fV0Filter->InitCutHistograms();
	fHistograms->Add(fV0Filter->GetCutHistograms());
  }
  

  SetUpCorrObjects();


  ///Set up ME histograms
  TList * MEHistograms = new TList();
  MEHistograms->SetName("MEHistograms");
  MEHistograms->SetOwner(kTRUE);
  fHistograms->Add(MEHistograms);

  hMETracks = new TObjArray();
  hMETracks->SetName("TrackArray");
  hMETracks->SetOwner(kTRUE);
  hMEPhotons = new TObjArray();
  hMEPhotons->SetName("PhotonArray");
  hMEPhotons->SetOwner(kTRUE);
  hMEPions = new TObjArray();
  hMEPions->SetName("PionArray");
  hMEPions->SetOwner(kTRUE);

  MEHistograms->Add(hMETracks);
  MEHistograms->Add(hMEPions);
  MEHistograms->Add(hMEPhotons);

  hMEvents = new TH2I("hMEvents", "Nevents vs centrality vertexz",
					  fAxisZ.GetNbins(), fAxisZ.GetBinLowEdge(1), fAxisZ.GetBinUpEdge(fAxisZ.GetNbins()),
 					  fAxisCent.GetNbins(), fAxisCent.GetBinLowEdge(1), fAxisCent.GetBinUpEdge(fAxisCent.GetNbins()));
  hMEvents->GetYaxis()->Set(fAxisCent.GetNbins(), fAxisCent.GetXbins()->GetArray());
  MEHistograms->Add(hMEvents);


  TList axesList;
  axesList.AddAt(&GetAxisEta(), 0);
  axesList.AddAt(&GetAxisPhi(), 1);
  axesList.AddAt(&GetAxisPt(), 2);
  axesList.SetOwner(kFALSE);
  
  TList piAxesList;
  piAxesList.AddAt(&GetAxisEta(), 0);
  piAxesList.AddAt(&GetAxisPhi(), 1);
  piAxesList.AddAt(&GetAxisPt(), 2);
  piAxesList.AddAt(&GetAxisPiMass(), 3);
  piAxesList.SetOwner(kFALSE);


  TList * outAxesList = new TList();
  outAxesList->Add(&fAxisCent);
  outAxesList->Add(&fAxisZ);

  for(Int_t iz = 0; iz < fAxisZ.GetNbins(); iz++) {
	TObjArray * trackArray = new TObjArray();
	trackArray->SetName(Form("METracks_%d", iz));
	trackArray->SetOwner(kTRUE);
	TObjArray * photonArray = new TObjArray();
	photonArray->SetName(Form("MEPhotons_%d", iz));
	photonArray->SetOwner(kTRUE);
	TObjArray * pionArray = new TObjArray();
	pionArray->SetName(Form("MEPions_%d", iz));
	pionArray->SetOwner(kTRUE);


	hMEPions->AddAt(pionArray, iz);
	hMETracks->AddAt(trackArray, iz);
	hMEPhotons->AddAt(photonArray, iz);

	for(Int_t ic = 0; ic < fAxisCent.GetNbins(); ic++) {

	  TString nameString = Form("%d_%d", iz, ic);
	  TString titleString = Form("%f < Z < %f ... %f cent %f", 
								 fAxisZ.GetBinLowEdge(iz+1), fAxisZ.GetBinUpEdge(iz+1), 
								 fAxisCent.GetBinLowEdge(ic+1), fAxisCent.GetBinUpEdge(ic+1));


	  THnSparseF * trackHistogram = CreateSparse(Form("tracks_%s", nameString.Data()), 
												 Form("tracks %s", titleString.Data()), &axesList );
	  trackArray->AddAt(trackHistogram, ic);

	  THnSparseF * photonHistogram = CreateSparse(Form("photons_%s", nameString.Data()), 
												 Form("photons %s", titleString.Data()), &axesList );
	  photonArray->AddAt(photonHistogram, ic);

	  THnSparseF * pionHistogram = CreateSparse(Form("pions_%s", nameString.Data()), 
												 Form("pions %s", titleString.Data()), &piAxesList );
	  pionArray->AddAt(pionHistogram, ic);
	}
  }

  PostData(1, fHistograms);
  PostData(2, fHistoGamma);
  PostData(3, fHistoPion);

}

///________________________________________________________________________
THnSparseF * AliAnalysisTaskdPhi::CreateSparse(TString nameString, TString titleString, TList * axesList) {
  ///Creat sparse
  const Int_t dim = axesList->GetSize();

  TAxis * axes[dim];
  Int_t bins[dim];
  Double_t min[dim];
  Double_t max[dim];

  for(Int_t i = 0; i<dim; i++) {
	TAxis * axis = dynamic_cast<TAxis*>(axesList->At(i));
	if(axis) {
	  axes[i] = axis;
  	} else {
	  cout << "AliAnalysisTaskdPhi::CreateSparse: Error error, all the axes are not present in axis list" << endl;
	  return NULL;
	}
  }

  for(Int_t i = 0; i<dim; i++) {
	bins[i] = axes[i]->GetNbins(); 
	min[i] = axes[i]->GetBinLowEdge(1);
	max[i] = axes[i]->GetBinUpEdge(axes[i]->GetNbins());
  }

  THnSparseF * sparse = new THnSparseF(Form("METracks_%s", nameString.Data()), 
											   Form("tracks %s", titleString.Data()), 
											   dim, bins, min, max);
  
  for(Int_t i = 0; i<dim; i++) {
	sparse->GetAxis(i)->SetNameTitle(axes[i]->GetName(), axes[i]->GetTitle() );
	if(axes[i]->GetXbins()->GetSize() > 0) {
	  sparse->SetBinEdges(i, axes[i]->GetXbins()->GetArray() );
	}
  }

  return sparse;
}

//________________________________________________________________________
void AliAnalysisTaskdPhi::UserExec(Option_t *) {
  ///User exec. 

  //if(! fV0Filter->EventIsSelected(fInputEvent)) return;


  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
  Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
  
  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
  if (!inputHandler) {
	cout << "cout no input event handler"<<endl;
	return;
  }

  
  if ( fV0Filter && !fV0Filter->GetPIDResponse() ) {
	if ( inputHandler->GetPIDResponse() ){
	  fV0Filter->SetPIDResponse( inputHandler->GetPIDResponse() );
	} else {
	  
	  //AOD case
	  if (isAOD){
		if (!fV0Filter->GetPIDResponse()){
		  fV0Filter->InitAODpidUtil(1);
		}
	  }
	}
  }

  Double_t centrality = 0.0;
  Double_t eventPlane = 0.0;
  Double_t vertexz = fInputEvent->GetPrimaryVertex()->GetZ();
  if(isAOD) {
    AliAODHeader * header = static_cast<AliAODHeader*>(fInputEvent->GetHeader());
	centrality = header->GetCentrality();
	eventPlane = header->GetEventplane();
  } else {
	centrality = static_cast<AliESDEvent*>(fInputEvent)->GetCentrality()->GetCentralityPercentile("kV0M");
	eventPlane = fInputEvent->GetEventplane()->GetEventplane("Q");
  }

  if(DebugLevel () > 15) {
	cout << "centrality: " << centrality <<  " " << GetBin(fAxisCent, centrality) << endl;
	cout << "vertexz: " << vertexz <<  " " << GetBin(fAxisZ, vertexz) << endl;
	cout << "eventPlane: " << eventPlane <<  " " << endl;
  }

  const Int_t centBin = GetBin(fAxisCent, centrality);
  const Int_t vertexBin = GetBin(fAxisZ, vertexz);

  if(centBin < 0 || vertexBin < 0) {
	AliError("bin out of range");
	return;
  }

  fGammas->Clear();
  fPions->Clear();

  TClonesArray * aodGammas = GetConversionGammas(isAOD);
  if(!aodGammas) {
	AliError("no aod gammas found!");
	return;
  }

  if(aodGammas->GetEntriesFast() > 0) {
	if( static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(0) == fL1 && 
		static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(1) == fL2 
		) {
	  return;
	}
	fL1 = static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(0);
	fL2 = static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(1);
	//cout << aodGammas->GetEntriesFast() << " " << fInputEvent->GetNumberOfTracks() << "c" << endl;
  }

  if(DebugLevel() > 1) printf("Number of conversion gammas %d \n", aodGammas->GetEntriesFast());
  for(Int_t ig = 0; ig < aodGammas->GetEntriesFast(); ig++) {
    AliAODConversionPhoton * photon = dynamic_cast<AliAODConversionPhoton*>(aodGammas->At(ig));
    
    if(!photon) continue;
    if(!fV0Filter || fV0Filter->PhotonIsSelected(static_cast<AliConversionPhotonBase*>(photon), fInputEvent)) {
      fGammas->Add(static_cast<TObject*>(photon));
    }
  }
  
  if(DebugLevel() > 4) printf("Number of accepted gammas %d \n", fGammas->GetEntriesFast());

  THnSparseF * trackMehist = GetMEHistogram(vertexBin, centBin, hMETracks); 
  hMEvents->Fill(vertexz, centrality);
  
  ///Add tracks to array
  TObjArray tracks;
  for(Int_t iTrack = 0; iTrack < fInputEvent->GetNumberOfTracks(); iTrack++) {

	AliVTrack * track = static_cast<AliVTrack*>(fInputEvent->GetTrack(iTrack));
	if(track->Pt() < 0.5) continue;
	if(TMath::Abs(track->Eta()) > 0.8) continue;
	tracks.Add(track);

	Double_t hValues[3] = { track->Eta(), track->Phi(), track->Pt() };
	if(trackMehist) trackMehist->Fill(hValues);
  }
  
  Process(fGammas, &tracks, vertexBin, centBin);

  PostData(1, fHistograms);
  PostData(2, fHistoGamma);
  PostData(3, fHistoPion);
  
}


//________________________________________________________________________
void AliAnalysisTaskdPhi::Process(TObjArray * gammas, TObjArray * tracks, Int_t vertexBin, Int_t centBin) {
  ///Process stuff

  if(DebugLevel() > 4) printf("Number of accepted tracks %d \n", tracks->GetEntriesFast());
 
  THnSparseF * gHisto = GetMEHistogram(vertexBin, centBin, hMEPhotons);
  THnSparseF * piHisto = GetMEHistogram(vertexBin, centBin, hMEPions);

  AliAnaConvCorrBase * gCorr = GetCorrObject(vertexBin, centBin, fPhotonCorr);
  AliAnaConvCorrBase * piCorr = GetCorrObject(vertexBin, centBin, fPionCorr);
  
  if(!gCorr || !piCorr) {
	AliError("corr object missing");
	return;
  }


  for(Int_t i1 = 0; i1 < gammas->GetEntriesFast(); i1++) {
   	AliAODConversionPhoton * ph1 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i1));
	Int_t tIDs[4] = {ph1->GetLabel(0), ph1->GetLabel(1), -1, -1};

	// AliESDv0 * v0 = ((AliESDEvent*) (fInputEvent))->GetV0(ph1->GetV0Index());
	// if(v0) {
	//   cout << v0->GetNindex() << " " << v0->GetPindex() << endl;
	//   cout << tIDs[1] << " " << tIDs[0] << endl;
	// }

	
	Bool_t leading = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs);

	if(leading) {
	  gCorr->CorrelateWithTracks( static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs, kFALSE);
	  Double_t phval[3] = {ph1->Eta(), ph1->Phi(), ph1->Pt()};
	  gHisto->Fill(phval);
	}

	for(Int_t i2 = 0; i2 < i1; i2++) {
	  AliAODConversionPhoton * ph2 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i2));

	  if( ph2->GetTrackLabelPositive()==ph1->GetTrackLabelPositive() 
		  || ph2->GetTrackLabelNegative()==ph1->GetTrackLabelNegative()
		  || ph2->GetTrackLabelNegative()==ph1->GetTrackLabelPositive()
		  || ph2->GetTrackLabelPositive()==ph1->GetTrackLabelNegative()) {
		continue;
	  }

	  AliAODConversionMother * pion = new AliAODConversionMother(ph1, ph2);
	  pion->SetLabels(i1, i2);

	  if(!fV0Filter || fV0Filter->MesonIsSelected(pion, kTRUE) ) {
		
		Bool_t leadingpi = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(pion), tracks, tIDs);
	
		if(leadingpi) {
		  //piCorr->FillTriggerCounters(pion);

		  Double_t pival[4] = { pion->Eta(), pion->Phi(), pion->Pt(), pion->M() };
		  piHisto->Fill(pival);
		
		  if(pion->M() > fAxisPiM.GetXmin() && pion->M() < fAxisPiM.GetXmax()) {
			tIDs[2] = ph2->GetLabel(0);
			tIDs[3] = ph2->GetLabel(1);
			piCorr->CorrelateWithTracks(pion, tracks, tIDs, kFALSE);
		  }
		} 
	  }
	}
  }
}


//________________________________________________________________________
void AliAnalysisTaskdPhi::Terminate(Option_t *) {
 
  // Draw result to the screen
  // Called once at the end of the query
}

//________________________________________________________________________
TClonesArray * AliAnalysisTaskdPhi::GetConversionGammas(Bool_t isAOD) {
  

  if(isAOD) {

	TClonesArray * gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
	if(gammas) {
	  return gammas;
	}
	//If not found try to locate branch
	FindDeltaAODBranchName(AODEvent());
  //gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
  //return gammas;
	
	gammas = dynamic_cast<TClonesArray*>(AODEvent()->FindListObject(fDeltaAODBranchName.Data()));
	return gammas;

  } else {
	TClonesArray * gammas = dynamic_cast<TClonesArray*>(GetInputData(1));
	return gammas;
  }

}

//________________________________________________________________________
void AliAnalysisTaskdPhi::FindDeltaAODBranchName(AliAODEvent * event){
  ///Find aod branch
  TList *list=event->GetList();
  for(Int_t ii=0;ii<list->GetEntries();ii++){
	TString name((list->At(ii))->GetName());
	if(name.BeginsWith("GammaConv")&&name.EndsWith("gamma")){
	  fDeltaAODBranchName=name;
	  AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
	  return;
	}
  }
}
Commit	Line	Data
3b77b2d1	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Authors: Svein Lindal *
	5	* Version 1.0 *
	6	* *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17	////////////////////////////////////////////////
	18	//---------------------------------------------
	19	// Class doing conversion gamma dPhi correlations
	20	// Gamma Conversion analysis
	21	//---------------------------------------------
	22	////////////////////////////////////////////////
	23
	24	#include "AliAnalysisTaskdPhi.h"
	25
	26	#include <TH2I.h>
	27	#include <TList.h>
	28	#include <TChain.h>
	29
	30	#include <AliAnalysisManager.h>
	31	#include <AliInputEventHandler.h>
	32	#include <AliESDInputHandler.h>
	33	#include <AliAODInputHandler.h>
	34	#include <AliAnalysisFilter.h>
	35
	36	#include "AliConversionCuts.h"
	37	#include "AliAODConversionPhoton.h"
	38	#include "AliAODConversionMother.h"
	39	#include "AliAnaConvCorrPhoton.h"
	40	#include "AliAnaConvCorrPion.h"
	41	#include "AliAnaConvIsolation.h"
	42	// Author Svein Lindal <slindal@fys.uio.no>
	43	using namespace std;
	44
	45	ClassImp(AliAnalysisTaskdPhi)
	46
	47
	48	//________________________________________________________________________
	49	AliAnalysisTaskdPhi::AliAnalysisTaskdPhi(const char *name) : AliAnalysisTaskSE(name),
	50	fHistograms(NULL),
	51	fHistoGamma(NULL),
	52	fHistoPion(NULL),
	53	fDielV0TrackFilter(NULL),
	54	fDielV0Filter(NULL),
	55	fDielPi0Filter(NULL),
	56	fDielTrackFilter(NULL),
	57	fV0Filter(NULL),
	58	fGammas(NULL),
	59	fPions(NULL),
	60	hMETracks(NULL),
	61	hMEPhotons(NULL),
	62	hMEPions(NULL),
	63	hMEvents(NULL),
	64	fPhotonCorr(NULL),
65	fPionCorr(NULL),
66	fIsoAna(NULL),
67	fL1(-1),
68	fL2(-1),
69	fDeltaAODBranchName("AliAODGammaConversion_gamma"),
70	fAxisPt(),
71	fAxisEta(),
72	fAxisPhi(),
73	fAxisCent(),
74	fAxisZ(),
75	fAxisPiM()
76	{
77	//constructor
78
79
80	fAxisPt.SetNameTitle("PtAxis", "Pt");
81	fAxisPt.Set(20, 0, 100);
82	fAxisEta.SetNameTitle("EtaAxis", "Eta");
83	fAxisEta.Set(160, -0.8, 0.8);
84	fAxisPhi.SetNameTitle("PhiAxis", "Phi");
85	fAxisPhi.Set(128, 0, TMath::TwoPi());
86	fAxisZ.SetNameTitle("ZAxis", "Z");
87	fAxisZ.Set(4, -10, 10);
88	fAxisCent.SetNameTitle("CentAxis", "Cent");
89
90	Double_t centbins[5] = {0, 10, 30, 60, 100.1};
91	fAxisCent.Set(4, centbins);
92
93
94	Double_t mbins[7] = {0.1, 0.11, 0.12, 0.15, 0.16, 0.18, 0.2};
95	fAxisPiM.SetNameTitle("InvMassPi0", "Invariant mass");
96	fAxisPiM.Set(6, mbins);
97
98	fGammas = new TObjArray();
99	fGammas->SetOwner(kFALSE);
100
101	fPions = new TObjArray();
102	fPions->SetOwner(kFALSE);
103
104	// Define input and output slots here
105	DefineInput(0, TChain::Class());
106	DefineInput(1, TClonesArray::Class());
107	DefineOutput(1, TList::Class());
108	DefineOutput(2, TList::Class());
109	DefineOutput(3, TList::Class());
110	}
111
112
113
114	//________________________________________________________________________
115	AliAnalysisTaskdPhi::~AliAnalysisTaskdPhi(){
116	//destructor
117	if(fPions)
118	delete fPions;
119	fPions = NULL;
120
121	if(fGammas)
122	delete fGammas;
123	fGammas = NULL;
124
125	if(fIsoAna)
126	delete fIsoAna;
127	fIsoAna = NULL;
128
129	if(fV0Filter)
130	delete fV0Filter;
131	fV0Filter = NULL;
132
133	if(fHistograms)
134	delete fHistograms;
135	fHistograms = NULL;
136
137	if(fHistoPion)
138	delete fHistoPion;
139	fHistoPion = NULL;
140
141	if(fHistoGamma)
142	delete fHistoGamma;
143	fHistoGamma = NULL;
144
145	}
146
147	///________________________________________________________________________
148	void AliAnalysisTaskdPhi::SetUpCorrObjects() {
149	///Creat corr obj
150	fIsoAna = new AliAnaConvIsolation();
151
152
153	fPhotonCorr = new TObjArray();
154	fPionCorr = new TObjArray();
155
156	TList * hPhoton = new TList();
157	hPhoton->SetName("hPhotonCorr");
158	hPhoton->SetOwner(kTRUE);
159	fHistoGamma->Add(hPhoton);
160
161	TList * hPion = new TList();
162	hPion->SetName("hPionCorr");
163	hPion->SetOwner(kTRUE);
164	fHistoPion->Add(hPion);
165
166
167	for(Int_t iz = 0; iz < fAxisZ.GetNbins(); iz++) {
168	TObjArray * photonArray = new TObjArray();
169	photonArray->SetOwner(kTRUE);
170	fPhotonCorr->AddAt(photonArray, iz);
171
172	TObjArray * pionArray = new TObjArray();
173	pionArray->SetOwner(kTRUE);
174	fPionCorr->AddAt(pionArray, iz);
175
176	TList * photonList = new TList();
177	photonList->SetName(Form("photon_%d", iz));
178	photonList->SetOwner(kTRUE);
179	hPhoton->AddAt(photonList, iz);
180
181	TList * pionList = new TList();
182	pionList->SetName(Form("pion_%d", iz));
183	pionList->SetOwner(kTRUE);
184	hPion->AddAt(pionList, iz);
185
186	for(Int_t ic = 0; ic < fAxisCent.GetNbins(); ic++) {
187
188	TString nameString = Form("%d_%d", iz, ic);
189	TString titleString = Form("%f < Z < %f ... %f cent %f",
190	fAxisZ.GetBinLowEdge(iz+1), fAxisZ.GetBinUpEdge(iz+1),
191	fAxisCent.GetBinLowEdge(ic+1), fAxisCent.GetBinUpEdge(ic+1));
192
193
194
195	AliAnaConvCorrPhoton * photonCorr = new AliAnaConvCorrPhoton(Form("PhotonCorr_%s", nameString.Data()), Form("photon %s", titleString.Data()));
196	photonArray->AddAt(photonCorr, ic);
197	photonCorr->GetAxistPt().Set(fAxisPt.GetNbins(), fAxisPt.GetXbins()->GetArray());
198	photonCorr->GetAxiscPt().Set(fAxisPt.GetNbins(), fAxisPt.GetXbins()->GetArray());
199	photonCorr->CreateHistograms();
200	photonList->Add(photonCorr->GetHistograms());
201
202	AliAnaConvCorrPion * pionCorr = new AliAnaConvCorrPion(Form("PionCorr_%s", nameString.Data()), Form("pion %s", titleString.Data()));
203	pionArray->AddAt(pionCorr, ic);
204	pionCorr->GetAxistPt().Set(fAxisPt.GetNbins(), fAxisPt.GetXbins()->GetArray());
205	pionCorr->GetAxiscPt().Set(fAxisPt.GetNbins(), fAxisPt.GetXbins()->GetArray());
206	pionCorr->GetAxisM().Set(fAxisPiM.GetNbins(), fAxisPiM.GetXbins()->GetArray());
207	pionCorr->CreateHistograms();
208	pionList->Add(pionCorr->GetHistograms());
209	}
210	}
211	}
212
213	//________________________________________________________________________
214	void AliAnalysisTaskdPhi::UserCreateOutputObjects() {
215	// Create histograms
216
217	fHistograms = new TList();
218	fHistograms->SetName("dPhi_histograms");
219	fHistograms->SetOwner(kTRUE);
220
221	fHistoGamma = new TList();
222	fHistoGamma->SetName("Gamma_histo");
223	fHistoGamma->SetOwner(kTRUE);
224
225	fHistoPion = new TList();
226	fHistoPion->SetName("Pion_histo");
227	fHistoPion->SetOwner(kTRUE);
228
229
230	if(fV0Filter) {
231	fV0Filter->InitCutHistograms();
232	fHistograms->Add(fV0Filter->GetCutHistograms());
233	}
234
235
236	SetUpCorrObjects();
237
238
239	///Set up ME histograms
240	TList * MEHistograms = new TList();
241	MEHistograms->SetName("MEHistograms");
242	MEHistograms->SetOwner(kTRUE);
243	fHistograms->Add(MEHistograms);
244
245	hMETracks = new TObjArray();
246	hMETracks->SetName("TrackArray");
247	hMETracks->SetOwner(kTRUE);
248	hMEPhotons = new TObjArray();
249	hMEPhotons->SetName("PhotonArray");
250	hMEPhotons->SetOwner(kTRUE);
251	hMEPions = new TObjArray();
252	hMEPions->SetName("PionArray");
253	hMEPions->SetOwner(kTRUE);
254
255	MEHistograms->Add(hMETracks);
256	MEHistograms->Add(hMEPions);
257	MEHistograms->Add(hMEPhotons);
258
259	hMEvents = new TH2I("hMEvents", "Nevents vs centrality vertexz",
260	fAxisZ.GetNbins(), fAxisZ.GetBinLowEdge(1), fAxisZ.GetBinUpEdge(fAxisZ.GetNbins()),
261	fAxisCent.GetNbins(), fAxisCent.GetBinLowEdge(1), fAxisCent.GetBinUpEdge(fAxisCent.GetNbins()));
262	hMEvents->GetYaxis()->Set(fAxisCent.GetNbins(), fAxisCent.GetXbins()->GetArray());
263	MEHistograms->Add(hMEvents);
264
265
266	TList axesList;
267	axesList.AddAt(&GetAxisEta(), 0);
268	axesList.AddAt(&GetAxisPhi(), 1);
269	axesList.AddAt(&GetAxisPt(), 2);
270	axesList.SetOwner(kFALSE);
271
272	TList piAxesList;
273	piAxesList.AddAt(&GetAxisEta(), 0);
274	piAxesList.AddAt(&GetAxisPhi(), 1);
275	piAxesList.AddAt(&GetAxisPt(), 2);
276	piAxesList.AddAt(&GetAxisPiMass(), 3);
277	piAxesList.SetOwner(kFALSE);
278
279
280	TList * outAxesList = new TList();
281	outAxesList->Add(&fAxisCent);
282	outAxesList->Add(&fAxisZ);
283
284	for(Int_t iz = 0; iz < fAxisZ.GetNbins(); iz++) {
285	TObjArray * trackArray = new TObjArray();
286	trackArray->SetName(Form("METracks_%d", iz));
287	trackArray->SetOwner(kTRUE);
288	TObjArray * photonArray = new TObjArray();
289	photonArray->SetName(Form("MEPhotons_%d", iz));
290	photonArray->SetOwner(kTRUE);
291	TObjArray * pionArray = new TObjArray();
292	pionArray->SetName(Form("MEPions_%d", iz));
293	pionArray->SetOwner(kTRUE);
294
295
296	hMEPions->AddAt(pionArray, iz);
297	hMETracks->AddAt(trackArray, iz);
298	hMEPhotons->AddAt(photonArray, iz);
299
300	for(Int_t ic = 0; ic < fAxisCent.GetNbins(); ic++) {
301
302	TString nameString = Form("%d_%d", iz, ic);
303	TString titleString = Form("%f < Z < %f ... %f cent %f",
304	fAxisZ.GetBinLowEdge(iz+1), fAxisZ.GetBinUpEdge(iz+1),
305	fAxisCent.GetBinLowEdge(ic+1), fAxisCent.GetBinUpEdge(ic+1));
306
307
308	THnSparseF * trackHistogram = CreateSparse(Form("tracks_%s", nameString.Data()),
309	Form("tracks %s", titleString.Data()), &axesList );
310	trackArray->AddAt(trackHistogram, ic);
311
312	THnSparseF * photonHistogram = CreateSparse(Form("photons_%s", nameString.Data()),
313	Form("photons %s", titleString.Data()), &axesList );
314	photonArray->AddAt(photonHistogram, ic);
315
316	THnSparseF * pionHistogram = CreateSparse(Form("pions_%s", nameString.Data()),
317	Form("pions %s", titleString.Data()), &piAxesList );
318	pionArray->AddAt(pionHistogram, ic);
319	}
320	}
321
322	PostData(1, fHistograms);
323	PostData(2, fHistoGamma);
324	PostData(3, fHistoPion);
325
326	}
327
328	///________________________________________________________________________
329	THnSparseF * AliAnalysisTaskdPhi::CreateSparse(TString nameString, TString titleString, TList * axesList) {
330	///Creat sparse
331	const Int_t dim = axesList->GetSize();
332
333	TAxis * axes[dim];
334	Int_t bins[dim];
335	Double_t min[dim];
336	Double_t max[dim];
337
338	for(Int_t i = 0; i<dim; i++) {
339	TAxis * axis = dynamic_cast<TAxis*>(axesList->At(i));
340	if(axis) {
341	axes[i] = axis;
342	} else {
343	cout << "AliAnalysisTaskdPhi::CreateSparse: Error error, all the axes are not present in axis list" << endl;
344	return NULL;
345	}
346	}
347
348	for(Int_t i = 0; i<dim; i++) {
349	bins[i] = axes[i]->GetNbins();
350	min[i] = axes[i]->GetBinLowEdge(1);
351	max[i] = axes[i]->GetBinUpEdge(axes[i]->GetNbins());
352	}
353
354	THnSparseF * sparse = new THnSparseF(Form("METracks_%s", nameString.Data()),
355	Form("tracks %s", titleString.Data()),
356	dim, bins, min, max);
357
358	for(Int_t i = 0; i<dim; i++) {
359	sparse->GetAxis(i)->SetNameTitle(axes[i]->GetName(), axes[i]->GetTitle() );
360	if(axes[i]->GetXbins()->GetSize() > 0) {
361	sparse->SetBinEdges(i, axes[i]->GetXbins()->GetArray() );
362	}
363	}
364
365	return sparse;
366	}
367
368	//________________________________________________________________________
369	void AliAnalysisTaskdPhi::UserExec(Option_t *) {
370	///User exec.
371
372	//if(! fV0Filter->EventIsSelected(fInputEvent)) return;
373
374
375	AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
376	Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
377
378	AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
379	if (!inputHandler) {
380	cout << "cout no input event handler"<<endl;
381	return;
382	}
383
384
385	if ( fV0Filter && !fV0Filter->GetPIDResponse() ) {
386	if ( inputHandler->GetPIDResponse() ){
387	fV0Filter->SetPIDResponse( inputHandler->GetPIDResponse() );
388	} else {
389
390	//AOD case
391	if (isAOD){
392	if (!fV0Filter->GetPIDResponse()){
393	fV0Filter->InitAODpidUtil(1);
394	}
395	}
396	}
397	}
398
399	Double_t centrality = 0.0;
400	Double_t eventPlane = 0.0;
401	Double_t vertexz = fInputEvent->GetPrimaryVertex()->GetZ();
402	if(isAOD) {
403	AliAODHeader * header = static_cast<AliAODHeader*>(fInputEvent->GetHeader());
404	centrality = header->GetCentrality();
405	eventPlane = header->GetEventplane();
406	} else {
407	centrality = static_cast<AliESDEvent*>(fInputEvent)->GetCentrality()->GetCentralityPercentile("kV0M");
408	eventPlane = fInputEvent->GetEventplane()->GetEventplane("Q");
409	}
410
411	if(DebugLevel () > 15) {
412	cout << "centrality: " << centrality << " " << GetBin(fAxisCent, centrality) << endl;
413	cout << "vertexz: " << vertexz << " " << GetBin(fAxisZ, vertexz) << endl;
414	cout << "eventPlane: " << eventPlane << " " << endl;
415	}
416
417	const Int_t centBin = GetBin(fAxisCent, centrality);
418	const Int_t vertexBin = GetBin(fAxisZ, vertexz);
419
420	if(centBin < 0 \|\| vertexBin < 0) {
421	AliError("bin out of range");
422	return;
423	}
424
425	fGammas->Clear();
426	fPions->Clear();
427
428	TClonesArray * aodGammas = GetConversionGammas(isAOD);
429	if(!aodGammas) {
430	AliError("no aod gammas found!");
431	return;
432	}
433
434	if(aodGammas->GetEntriesFast() > 0) {
435	if( static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(0) == fL1 &&
436	static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(1) == fL2
437	) {
438	return;
439	}
440	fL1 = static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(0);
441	fL2 = static_cast<AliAODConversionParticle*>(aodGammas->At(0))->GetLabel(1);
442	//cout << aodGammas->GetEntriesFast() << " " << fInputEvent->GetNumberOfTracks() << "c" << endl;
443	}
444
445	if(DebugLevel() > 1) printf("Number of conversion gammas %d \n", aodGammas->GetEntriesFast());
446	for(Int_t ig = 0; ig < aodGammas->GetEntriesFast(); ig++) {
447	AliAODConversionPhoton * photon = dynamic_cast<AliAODConversionPhoton*>(aodGammas->At(ig));
448
449	if(!photon) continue;
450	if(!fV0Filter \|\| fV0Filter->PhotonIsSelected(static_cast<AliConversionPhotonBase*>(photon), fInputEvent)) {
451	fGammas->Add(static_cast<TObject*>(photon));
452	}
453	}
454
455	if(DebugLevel() > 4) printf("Number of accepted gammas %d \n", fGammas->GetEntriesFast());
456
457	THnSparseF * trackMehist = GetMEHistogram(vertexBin, centBin, hMETracks);
458	hMEvents->Fill(vertexz, centrality);
459
460	///Add tracks to array
461	TObjArray tracks;
462	for(Int_t iTrack = 0; iTrack < fInputEvent->GetNumberOfTracks(); iTrack++) {
463
464	AliVTrack * track = static_cast<AliVTrack*>(fInputEvent->GetTrack(iTrack));
465	if(track->Pt() < 0.5) continue;
466	if(TMath::Abs(track->Eta()) > 0.8) continue;
467	tracks.Add(track);
468
469	Double_t hValues[3] = { track->Eta(), track->Phi(), track->Pt() };
470	if(trackMehist) trackMehist->Fill(hValues);
471	}
472
473	Process(fGammas, &tracks, vertexBin, centBin);
474
475	PostData(1, fHistograms);
476	PostData(2, fHistoGamma);
477	PostData(3, fHistoPion);
478
479	}
480
481
482	//________________________________________________________________________
483	void AliAnalysisTaskdPhi::Process(TObjArray * gammas, TObjArray * tracks, Int_t vertexBin, Int_t centBin) {
484	///Process stuff
485
486	if(DebugLevel() > 4) printf("Number of accepted tracks %d \n", tracks->GetEntriesFast());
487
488	THnSparseF * gHisto = GetMEHistogram(vertexBin, centBin, hMEPhotons);
489	THnSparseF * piHisto = GetMEHistogram(vertexBin, centBin, hMEPions);
490
491	AliAnaConvCorrBase * gCorr = GetCorrObject(vertexBin, centBin, fPhotonCorr);
492	AliAnaConvCorrBase * piCorr = GetCorrObject(vertexBin, centBin, fPionCorr);
493
494	if(!gCorr \|\| !piCorr) {
495	AliError("corr object missing");
496	return;
497	}
498
499
500
501	for(Int_t i1 = 0; i1 < gammas->GetEntriesFast(); i1++) {
502	AliAODConversionPhoton * ph1 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i1));
503	Int_t tIDs[4] = {ph1->GetLabel(0), ph1->GetLabel(1), -1, -1};
504
505	// AliESDv0 * v0 = ((AliESDEvent*) (fInputEvent))->GetV0(ph1->GetV0Index());
506	// if(v0) {
507	// cout << v0->GetNindex() << " " << v0->GetPindex() << endl;
508	// cout << tIDs[1] << " " << tIDs[0] << endl;
509	// }
510
511
512	Bool_t leading = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs);
513
514	if(leading) {
515	gCorr->CorrelateWithTracks( static_cast<AliAODConversionParticle*>(ph1), tracks, tIDs, kFALSE);
516	Double_t phval[3] = {ph1->Eta(), ph1->Phi(), ph1->Pt()};
517	gHisto->Fill(phval);
518	}
519
520	for(Int_t i2 = 0; i2 < i1; i2++) {
521	AliAODConversionPhoton * ph2 = static_cast<AliAODConversionPhoton*>(gammas->UncheckedAt(i2));
522
523	if( ph2->GetTrackLabelPositive()==ph1->GetTrackLabelPositive()
524	\|\| ph2->GetTrackLabelNegative()==ph1->GetTrackLabelNegative()
525	\|\| ph2->GetTrackLabelNegative()==ph1->GetTrackLabelPositive()
526	\|\| ph2->GetTrackLabelPositive()==ph1->GetTrackLabelNegative()) {
527	continue;
528	}
529
530	AliAODConversionMother * pion = new AliAODConversionMother(ph1, ph2);
531	pion->SetLabels(i1, i2);
532
533	if(!fV0Filter \|\| fV0Filter->MesonIsSelected(pion, kTRUE) ) {
534
535	Bool_t leadingpi = fIsoAna->IsLeading(static_cast<AliAODConversionParticle*>(pion), tracks, tIDs);
536
537	if(leadingpi) {
538	//piCorr->FillTriggerCounters(pion);
539
540	Double_t pival[4] = { pion->Eta(), pion->Phi(), pion->Pt(), pion->M() };
541	piHisto->Fill(pival);
542
543	if(pion->M() > fAxisPiM.GetXmin() && pion->M() < fAxisPiM.GetXmax()) {
544	tIDs[2] = ph2->GetLabel(0);
545	tIDs[3] = ph2->GetLabel(1);
546	piCorr->CorrelateWithTracks(pion, tracks, tIDs, kFALSE);
547	}
548	}
549	}
550	}
551	}
552	}
553
554
555	//________________________________________________________________________
556	void AliAnalysisTaskdPhi::Terminate(Option_t *) {
557
558	// Draw result to the screen
559	// Called once at the end of the query
560	}
561
562	//________________________________________________________________________
563	TClonesArray * AliAnalysisTaskdPhi::GetConversionGammas(Bool_t isAOD) {
564
565
566	if(isAOD) {
567
568	TClonesArray * gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
569	if(gammas) {
570	return gammas;
571	}
572	//If not found try to locate branch
573	FindDeltaAODBranchName(AODEvent());
574	//gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
575	//return gammas;
576
577	gammas = dynamic_cast<TClonesArray*>(AODEvent()->FindListObject(fDeltaAODBranchName.Data()));
578	return gammas;
579
580	} else {
581	TClonesArray * gammas = dynamic_cast<TClonesArray*>(GetInputData(1));
582	return gammas;
583	}
584
585	}
586
587	//________________________________________________________________________
588	void AliAnalysisTaskdPhi::FindDeltaAODBranchName(AliAODEvent * event){
589	///Find aod branch
590	TList *list=event->GetList();
591	for(Int_t ii=0;ii<list->GetEntries();ii++){
592	TString name((list->At(ii))->GetName());
593	if(name.BeginsWith("GammaConv")&&name.EndsWith("gamma")){
594	fDeltaAODBranchName=name;
595	AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
596	return;
597	}
598	}
599	}
600
601