[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliHFCorrelator.cxx

/**************************************************************************\r
 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
 *                                                                        *\r
 * Author: The ALICE Off-line Project.                                    *\r
 * Contributors are mentioned in the code where appropriate.              *\r
 *                                                                        *\r
 * Permission to use, copy, modify and distribute this software and its   *\r
 * documentation strictly for non-commercial purposes is hereby granted   *\r
 * without fee, provided that the above copyright notice appears in all   *\r
 * copies and that both the copyright notice and this permission notice   *\r
 * appear in the supporting documentation. The authors make no claims     *\r
 * about the suitability of this software for any purpose. It is          *\r
 * provided "as is" without express or implied warranty.                  *\r
 **************************************************************************/\r
//\r
//\r
//             Base class for Heavy Flavour Correlations Analysis\r
//             Single Event and Mixed Event Analysis are implemented\r
//\r
//-----------------------------------------------------------------------\r
//          \r
//\r
//						   Author S.Bjelogrlic\r
//                         Utrecht University \r
//                      sandro.bjelogrlic@cern.ch\r
//\r
//-----------------------------------------------------------------------\r
\r
/* $Id$ */\r
\r
#include <TParticle.h>\r
#include <TVector3.h>\r
#include <TChain.h>\r
#include "TROOT.h"\r
#include "AliHFCorrelator.h"\r
#include "AliRDHFCutsDStartoKpipi.h"\r
#include "AliHFAssociatedTrackCuts.h"\r
#include "AliEventPoolManager.h"\r
#include "AliReducedParticle.h"\r
#include "AliCentrality.h"\r
#include "AliAODMCParticle.h"\r
\r
using std::cout;\r
using std::endl;\r
\r
//_____________________________________________________\r
AliHFCorrelator::AliHFCorrelator() :\r
//\r
// default constructor\r
//\r
TNamed(),\r
fPoolMgr(0x0),         \r
fPool(0x0),\r
fhadcuts(0x0),\r
fAODEvent(0x0),\r
fAssociatedTracks(0x0),\r
fmcArray(0x0),\r
fReducedPart(0x0),\r
fD0cand(0x0), \r
fhypD0(0), \r
fDCharge(0),\r
\r
fmixing(kFALSE),\r
fmontecarlo(kFALSE),\r
fsystem(kFALSE),\r
fUseReco(kTRUE),\r
fselect(kUndefined),\r
\r
fUseImpactParameter(0),\r
fPIDmode(0),\r
\r
fNofTracks(0),\r
fPoolContent(0),\r
\r
fPhiMin(0),\r
fPhiMax(0),\r
\r
fPtTrigger(0),\r
fPhiTrigger(0),\r
fEtaTrigger(0),\r
\r
\r
fDeltaPhi(0),\r
fDeltaEta(0),\r
fk0InvMass(0)\r
\r
{\r
	// default constructor	\r
}\r
\r
\r
\r
//_____________________________________________________\r
AliHFCorrelator::AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t ppOrPbPb) :\r
TNamed(name,"title"),\r
fPoolMgr(0x0),         \r
fPool(0x0),\r
fhadcuts(0x0),\r
fAODEvent(0x0),\r
fAssociatedTracks(0x0),\r
fmcArray(0x0),\r
fReducedPart(0x0),\r
fD0cand(0x0), \r
fhypD0(0),\r
fDCharge(0),\r
\r
fmixing(kFALSE),\r
fmontecarlo(kFALSE),\r
fsystem(ppOrPbPb),\r
fUseReco(kTRUE),\r
fselect(kUndefined),\r
fUseImpactParameter(0),\r
fPIDmode(0),\r
\r
fNofTracks(0),\r
fPoolContent(0),\r
\r
fPhiMin(0),\r
fPhiMax(0),\r
\r
fPtTrigger(0),\r
fPhiTrigger(0),\r
fEtaTrigger(0),\r
\r
\r
fDeltaPhi(0),\r
fDeltaEta(0),\r
fk0InvMass(0)\r
{\r
	fhadcuts = cuts; \r
}\r
\r
//_____________________________________________________\r
AliHFCorrelator::~AliHFCorrelator() \r
{\r
//\r
// destructor\r
//	\r
	\r
	if(fPoolMgr)  {delete fPoolMgr; fPoolMgr=0;}       \r
	if(fPool) {delete fPool; fPool=0;}\r
	if(fhadcuts) {delete fhadcuts; fhadcuts=0;}\r
	if(fAODEvent) {delete fAODEvent; fAODEvent=0;}\r
	if(fAssociatedTracks) {delete fAssociatedTracks; fAssociatedTracks=0;}\r
	if(fmcArray) {delete fmcArray; fmcArray=0;}\r
	if(fReducedPart) {delete fReducedPart; fReducedPart=0;}\r
	if(fD0cand) {delete fD0cand; fD0cand=0;}\r
	\r
	\r
	if(fNofTracks) fNofTracks = 0;\r
	\r
	if(fPhiMin) fPhiMin = 0;\r
	if(fPhiMax) fPhiMax = 0;\r
	\r
	if(fPtTrigger) fPtTrigger=0;\r
	if(fPhiTrigger) fPhiTrigger=0;\r
	if(fEtaTrigger) fEtaTrigger=0;\r
	\r
	if(fDeltaPhi) fDeltaPhi=0;\r
	if(fDeltaEta) fDeltaEta=0;\r
	\r
	if(fk0InvMass) fk0InvMass=0;\r
	\r
}\r
//_____________________________________________________\r
Bool_t AliHFCorrelator::DefineEventPool(){\r
	// definition of the Pool Manager for Event Mixing\r
	\r
\r
	Int_t MaxNofEvents = fhadcuts->GetMaxNEventsInPool();\r
	Int_t MinNofTracks = fhadcuts->GetMinNTracksInPool();\r
	Int_t NofCentBins = fhadcuts->GetNCentPoolBins();\r
	Double_t * CentBins = fhadcuts->GetCentPoolBins();\r
	Int_t NofZVrtxBins = fhadcuts->GetNZvtxPoolBins();\r
	Double_t *ZVrtxBins = fhadcuts->GetZvtxPoolBins();\r
		\r
			\r
	fPoolMgr = new AliEventPoolManager(MaxNofEvents, MinNofTracks, NofCentBins, CentBins, NofZVrtxBins, ZVrtxBins);\r
	if(!fPoolMgr) return kFALSE;\r
	return kTRUE;\r
}\r
//_____________________________________________________\r
Bool_t AliHFCorrelator::Initialize(){\r
	\r
    //  std::cout << "AliHFCorrelator::Initialize"<< std::endl;\r
//  AliInfo("AliHFCorrelator::Initialize") ;\r
  if(!fAODEvent){\r
    AliInfo("No AOD event") ;\r
    return kFALSE;\r
  }\r
    //std::cout << "No AOD event" << std::endl;\r
	\r
	AliCentrality *centralityObj = 0;\r
	//Int_t multiplicity = -1;\r
	Double_t MultipOrCent = -1;\r
	\r
	// initialize the pool for event mixing\r
	if(!fsystem){ // pp\r
	//multiplicity = fAODEvent->GetNTracks();\r
        MultipOrCent = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(fAODEvent,-1.,1.);\r
	//	MultipOrCent = multiplicity; // convert from Int_t to Double_t\r
     //   AliInfo(Form("Multiplicity is %f", MultipOrCent));\r
	}\r
	if(fsystem){ // PbPb\r
		\r
		centralityObj = fAODEvent->GetHeader()->GetCentralityP();\r
		MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");\r
//		AliInfo(Form("Centrality is %f", MultipOrCent));\r
	}\r
	\r
	AliAODVertex *vtx = fAODEvent->GetPrimaryVertex();\r
	Double_t zvertex = vtx->GetZ(); // zvertex\r
	Double_t * CentBins = fhadcuts->GetCentPoolBins();\r
	Double_t poolmin=CentBins[0];\r
	Double_t poolmax=CentBins[fhadcuts->GetNCentPoolBins()];\r
\r
	\r
		if(TMath::Abs(zvertex)>=10 || MultipOrCent>poolmax || MultipOrCent < poolmin) {\r
		if(!fsystem)AliInfo(Form("pp or pA Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,MultipOrCent));\r
		if(fsystem) AliInfo(Form("PbPb Event with Zvertex = %.2f cm and centrality = %.1f  out of pool bounds, SKIPPING",zvertex,MultipOrCent));\r
\r
			return kFALSE;\r
		}\r
	\r
	fPool = fPoolMgr->GetEventPool(MultipOrCent, zvertex);\r
	\r
	if (!fPool){\r
		AliInfo(Form("No pool found for multiplicity = %f, zVtx = %f cm", MultipOrCent, zvertex));\r
	    return kFALSE;\r
	}\r
	//fPool->PrintInfo();\r
	return kTRUE;\r
}\r
\r
//_____________________________________________________\r
Bool_t AliHFCorrelator::ProcessEventPool(){\r
	 // analysis on Mixed Events\r
	//cout << "AliHFCorrelator::ProcessEventPool"<< endl;\r
		if(!fmixing) return kFALSE;\r
		if(!fPool->IsReady()) return kFALSE;\r
		if(fPool->GetCurrentNEvents()<fhadcuts->GetMinEventsToMix()) return kFALSE;\r
	//	fPool->PrintInfo();\r
		fPoolContent = fPool->GetCurrentNEvents();\r
		\r
		return kTRUE;\r
	\r
}\r
\r
//_____________________________________________________\r
Bool_t AliHFCorrelator::ProcessAssociatedTracks(Int_t EventLoopIndex, const TObjArray* associatedTracks){\r
  // associatedTracks is not deleted, it should be (if needed) deleted in the user task\r
  \r
  if(!fmixing){ // analysis on Single Event\r
    if(fAssociatedTracks){\r
      fAssociatedTracks->Delete();\r
      delete fAssociatedTracks;\r
    }      \r
    if(fselect==kHadron || fselect ==kKaon){\r
      fAssociatedTracks = AcceptAndReduceTracks(fAODEvent);\r
      fAssociatedTracks->SetOwner(kTRUE);\r
    }\r
    if(fselect==kKZero) {\r
      fAssociatedTracks = AcceptAndReduceKZero(fAODEvent);\r
      fAssociatedTracks->SetOwner(kTRUE);\r
    }	\r
    if(fselect==kElectron && associatedTracks) {\r
      fAssociatedTracks=new TObjArray(*associatedTracks);// Maybe better to call the copy constructor\r
      fAssociatedTracks->SetOwner(kFALSE);\r
    }\r
    \r
  }\r
  \r
  if(fmixing) { // analysis on Mixed Events\r
		\r
			\r
    fAssociatedTracks = fPool->GetEvent(EventLoopIndex);\r
				\r
    \r
    \r
    \r
  } // end if mixing\r
  \r
  if(!fAssociatedTracks) return kFALSE;\r
  \r
  fNofTracks = fAssociatedTracks->GetEntriesFast(); \r
  \r
  return kTRUE;\r
	\r
}\r
//_____________________________________________________\r
Bool_t AliHFCorrelator::Correlate(Int_t loopindex){\r
\r
	if(loopindex >= fNofTracks) return kFALSE;\r
	if(!fAssociatedTracks) return kFALSE;\r
	\r
	fReducedPart = (AliReducedParticle*)fAssociatedTracks->At(loopindex);\r
	\r
\r
	fDeltaPhi = SetCorrectPhiRange(fPhiTrigger - fReducedPart->Phi());\r
	\r
	fDeltaEta = fEtaTrigger - fReducedPart->Eta();\r
\r
	return kTRUE;\r
	\r
}\r
		\r
//_____________________________________________________\r
Bool_t AliHFCorrelator::PoolUpdate(const TObjArray* associatedTracks){\r
\r
	if(!fmixing) return kFALSE;\r
	if(!fPool) return kFALSE;\r
	if(fmixing) { // update the pool for Event Mixing\r
		TObjArray* objArr = NULL;\r
		if(fselect==kHadron || fselect==kKaon) objArr = (TObjArray*)AcceptAndReduceTracks(fAODEvent);\r
		else if(fselect==kKZero) objArr = (TObjArray*)AcceptAndReduceKZero(fAODEvent);\r
		else if(fselect==kElectron && associatedTracks){\r
		  objArr = new TObjArray(*associatedTracks);\r
		}\r
		else return kFALSE;\r
		if(objArr->GetEntriesFast()>0) fPool->UpdatePool(objArr); // updating the pool only if there are entries in the array\r
	}\r
		\r
	return kTRUE;\r
	\r
}\r
		\r
//_____________________________________________________\r
Double_t AliHFCorrelator::SetCorrectPhiRange(Double_t phi){\r
	Double_t pi = TMath::Pi();\r
	\r
	if(phi<fPhiMin) phi = phi + 2*pi;\r
	if(phi>fPhiMax) phi = phi - 2*pi;\r
	\r
	return phi;\r
}\r
\r
//_____________________________________________________\r
TObjArray*  AliHFCorrelator::AcceptAndReduceTracks(AliAODEvent* inputEvent){\r
\r
  Double_t weight=1.;\r
  Int_t nTracks = inputEvent->GetNTracks();\r
  AliAODVertex * vtx = inputEvent->GetPrimaryVertex();\r
  Double_t pos[3],cov[6];\r
  vtx->GetXYZ(pos);\r
  vtx->GetCovarianceMatrix(cov);\r
  const AliESDVertex vESD(pos,cov,100.,100);\r
  \r
  Double_t Bz = inputEvent->GetMagneticField();\r
	\r
  \r
  TObjArray* tracksClone = new TObjArray;\r
  tracksClone->SetOwner(kTRUE);\r
  \r
  //*******************************************************\r
  // use reconstruction\r
  if(fUseReco){\r
    for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {\r
      AliAODTrack* track = inputEvent->GetTrack(iTrack);\r
      if (!track) continue;\r
      if(!fhadcuts->IsHadronSelected(track,&vESD,Bz)) continue; // apply ESD level selections\r
      if(!fhadcuts->Charge(fDCharge,track)) continue; // apply selection on charge, if required\r
\r
      Double_t pT = track->Pt();\r
      \r
      //compute impact parameter\r
      Double_t d0z0[2],covd0z0[3];\r
      Double_t d0=-999999.;\r
      if(fUseImpactParameter) track->PropagateToDCA(vtx,Bz,100,d0z0,covd0z0);\r
      else d0z0[0] = 1. ; // random number - be careful with the cuts you applied\r
      \r
      if(fUseImpactParameter==1) d0 = TMath::Abs(d0z0[0]); // use impact parameter\r
      if(fUseImpactParameter==2) { // use impact parameter over resolution\r
	if(TMath::Abs(covd0z0[0])>0.00000001) d0 = TMath::Abs(d0z0[0])/TMath::Sqrt(covd0z0[0]); \r
	else d0 = -1.; // if the resoultion is Zero, rejects the track - to be on the safe side\r
	\r
      }\r
      \r
      if(fmontecarlo) {// THIS TO BE CHECKED\r
        Int_t hadLabel = track->GetLabel();\r
	if(hadLabel < 0) continue;	\r
      }\r
      \r
      if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts\r
      Bool_t rejectsoftpi = kTRUE;// TO BE CHECKED: DO WE WANT IT TO kTRUE AS A DEFAULT?\r
      if(fD0cand && !fmixing) rejectsoftpi = fhadcuts->InvMassDstarRejection(fD0cand,track,fhypD0); // TO BE CHECKED: WHY NOT FOR EM?\r
      \r
      \r
      if(fselect ==kKaon){	\r
	if(!fhadcuts->CheckKaonCompatibility(track,fmontecarlo,fmcArray,fPIDmode)) continue; // check if it is a Kaon - data and MC\r
      }\r
      weight=fhadcuts->GetTrackWeight(pT,track->Eta(),pos[2]);\r
      tracksClone->Add(new AliReducedParticle(track->Eta(), track->Phi(), pT,track->GetLabel(),track->GetID(),d0,rejectsoftpi,track->Charge(),weight));\r
    } // end loop on tracks\r
  } // end if use reconstruction kTRUE\r
  \r
  //*******************************************************\r
  \r
  //use MC truth\r
  if(fmontecarlo && !fUseReco){\r
    \r
    for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) { \r
      AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));\r
      if (!mcPart) {\r
	AliWarning("MC Particle not found in tree, skipping"); \r
	continue;\r
      }\r
      if(!mcPart->Charge()) continue; // consider only charged tracks\r
      \r
      Int_t PDG =TMath::Abs(mcPart->PdgCode()); \r
if(fselect ==kHadron) {if(!((PDG==321)||(PDG==211)||(PDG==2212)||(PDG==13)||(PDG==11))) continue;} // select only if kaon, pion, proton, muon or electron\r
      else if(fselect ==kKaon) {if(!(PDG==321)) continue;} // select only if kaon\r
      else if(fselect ==kElectron) {if(!(PDG==11)) continue;} // select only if electron\r
\r
      Double_t pT = mcPart->Pt();\r
      Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet\r
      if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts\r
      \r
      tracksClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,iPart,-1,d0,kFALSE,mcPart->Charge()));\r
    }\r
    \r
  } // end if use  MC truth\r
  \r
  \r
  return tracksClone;\r
}\r
\r
//_____________________________________________________\r
TObjArray*  AliHFCorrelator::AcceptAndReduceKZero(AliAODEvent* inputEvent){\r
	\r
	Int_t nOfVZeros = inputEvent->GetNumberOfV0s();\r
	TObjArray* KZeroClone = new TObjArray;\r
	AliAODVertex *vertex1 = (AliAODVertex*)inputEvent->GetPrimaryVertex();\r
\r
 // use reconstruction	 	\r
     if(fUseReco){\r
	Int_t v0label = -1;\r
	Int_t pdgDgK0toPipi[2] = {211,211};\r
	Double_t mPDGK0=0.497614;//TDatabasePDG::Instance()->GetParticle(310)->Mass();\r
	const Int_t size = inputEvent->GetNumberOfV0s();\r
	Int_t prevnegID[size];\r
	Int_t prevposID[size];\r
	for(Int_t iv0 =0; iv0< nOfVZeros; iv0++){// loop on all v0 candidates\r
		AliAODv0 *v0 = (static_cast<AliAODEvent*> (inputEvent))->GetV0(iv0);\r
		if(!v0) {\r
		  AliInfo(Form("is not a v0 at step, %d", iv0)) ;\r
		  //cout << "is not a v0 at step " << iv0 << endl;\r
		  continue;\r
		}\r
		\r
		if(!fhadcuts->IsKZeroSelected(v0,vertex1)) continue; // check if it is a k0\r
	    \r
		// checks to avoid double counting\r
		Int_t negID = -999;\r
		Int_t posID = -998;\r
		//Int_t a = 0;\r
		prevnegID[iv0] = -997;\r
		prevposID[iv0]= -996;\r
		negID = v0->GetNegID();\r
		posID = v0->GetPosID();\r
		Bool_t DoubleCounting = kFALSE;\r
		\r
		for(Int_t k=0; k<iv0; k++){\r
			if(((negID==prevnegID[k])&&(posID==prevposID[k]))||((negID==prevposID[k])&&(posID==prevnegID[k]))){\r
				DoubleCounting = kTRUE;\r
				//a=k;\r
				break;\r
			}//end if\r
		} // end for\r
		\r
		if(DoubleCounting) continue;\r
		else{ \r
			prevposID[iv0] = posID; \r
			prevnegID[iv0] = negID;\r
		}\r
		\r
		if(fmontecarlo)	v0label = v0->MatchToMC(310,fmcArray, 0, pdgDgK0toPipi); //match a K0 short\r
		Double_t k0pt = v0->Pt();\r
		Double_t k0eta = v0->Eta();\r
		Double_t k0Phi = v0->Phi();\r
	    fk0InvMass = v0->MassK0Short();	\r
		\r
		//if (loopindex == 0) {\r
		//	if(!plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectra"))->Fill(k0InvMass,k0pt); // spectra for all k0\r
		//	if(plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectraifHF"))->Fill(k0InvMass,k0pt); // spectra for k0 in association with a D*\r
		//}\r
		// if there are more D* candidates per event, loopindex == 0 makes sure you fill the mass spectra only once!\r
		\r
		if(TMath::Abs(fk0InvMass-mPDGK0)>3*0.004) continue; // select candidates within 3 sigma\r
		KZeroClone->Add(new AliReducedParticle(k0eta,k0Phi,k0pt,v0label));\r
		\r
	}\r
     } // end if use reconstruction kTRUE\r
	\r
\r
\r
//*********************************************************************//\r
     //use MC truth\r
     if(fmontecarlo && !fUseReco){\r
		\r
		for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) { \r
			AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));\r
			if (!mcPart) {\r
				AliWarning("MC Particle not found in tree, skipping"); \r
				continue;\r
			}\r
			\r
			Int_t PDG =TMath::Abs(mcPart->PdgCode()); \r
			if(!(PDG==310)) continue; // select only if k0 short\r
		\r
			Double_t pT = mcPart->Pt();\r
            Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet\r
			if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts\r
			\r
			KZeroClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,iPart,-1,d0,kTRUE,mcPart->Charge()));\r
		}\r
		\r
	} // end if use  MC truth\r
\r
\r
\r
	return KZeroClone;\r
}\r
Commit	Line	Data
e9b1e63f	1	/**************************************************************************\r
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
	3	* *\r
	4	* Author: The ALICE Off-line Project. *\r
	5	* Contributors are mentioned in the code where appropriate. *\r
	6	* *\r
	7	* Permission to use, copy, modify and distribute this software and its *\r
	8	* documentation strictly for non-commercial purposes is hereby granted *\r
	9	* without fee, provided that the above copyright notice appears in all *\r
	10	* copies and that both the copyright notice and this permission notice *\r
	11	* appear in the supporting documentation. The authors make no claims *\r
	12	* about the suitability of this software for any purpose. It is *\r
	13	* provided "as is" without express or implied warranty. *\r
	14	**************************************************************************/\r
	15	//\r
	16	//\r
	17	// Base class for Heavy Flavour Correlations Analysis\r
	18	// Single Event and Mixed Event Analysis are implemented\r
	19	//\r
	20	//-----------------------------------------------------------------------\r
	21	// \r
	22	//\r
	23	// Author S.Bjelogrlic\r
	24	// Utrecht University \r
	25	// sandro.bjelogrlic@cern.ch\r
	26	//\r
	27	//-----------------------------------------------------------------------\r
	28	\r
	29	/* $Id$ */\r
	30	\r
	31	#include <TParticle.h>\r
	32	#include <TVector3.h>\r
	33	#include <TChain.h>\r
	34	#include "TROOT.h"\r
	35	#include "AliHFCorrelator.h"\r
	36	#include "AliRDHFCutsDStartoKpipi.h"\r
	37	#include "AliHFAssociatedTrackCuts.h"\r
	38	#include "AliEventPoolManager.h"\r
	39	#include "AliReducedParticle.h"\r
	40	#include "AliCentrality.h"\r
	41	#include "AliAODMCParticle.h"\r
	42	\r
	43	using std::cout;\r
	44	using std::endl;\r
	45	\r
	46	//_____________________________________________________\r
	47	AliHFCorrelator::AliHFCorrelator() :\r
	48	//\r
	49	// default constructor\r
	50	//\r
	51	TNamed(),\r
	52	fPoolMgr(0x0), \r
	53	fPool(0x0),\r
	54	fhadcuts(0x0),\r
	55	fAODEvent(0x0),\r
	56	fAssociatedTracks(0x0),\r
	57	fmcArray(0x0),\r
	58	fReducedPart(0x0),\r
	59	fD0cand(0x0), \r
	60	fhypD0(0), \r
	61	fDCharge(0),\r
	62	\r
	63	fmixing(kFALSE),\r
	64	fmontecarlo(kFALSE),\r
65	fsystem(kFALSE),\r
66	fUseReco(kTRUE),\r
67	fselect(kUndefined),\r
68	\r
69	fUseImpactParameter(0),\r
70	fPIDmode(0),\r
71	\r
72	fNofTracks(0),\r
73	fPoolContent(0),\r
74	\r
75	fPhiMin(0),\r
76	fPhiMax(0),\r
77	\r
78	fPtTrigger(0),\r
79	fPhiTrigger(0),\r
80	fEtaTrigger(0),\r
81	\r
82	\r
83	fDeltaPhi(0),\r
84	fDeltaEta(0),\r
85	fk0InvMass(0)\r
86	\r
87	{\r
88	// default constructor \r
89	}\r
90	\r
91	\r
92	\r
93	//_____________________________________________________\r
94	AliHFCorrelator::AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t ppOrPbPb) :\r
95	TNamed(name,"title"),\r
96	fPoolMgr(0x0), \r
97	fPool(0x0),\r
98	fhadcuts(0x0),\r
99	fAODEvent(0x0),\r
100	fAssociatedTracks(0x0),\r
101	fmcArray(0x0),\r
102	fReducedPart(0x0),\r
103	fD0cand(0x0), \r
104	fhypD0(0),\r
105	fDCharge(0),\r
106	\r
107	fmixing(kFALSE),\r
108	fmontecarlo(kFALSE),\r
109	fsystem(ppOrPbPb),\r
110	fUseReco(kTRUE),\r
111	fselect(kUndefined),\r
112	fUseImpactParameter(0),\r
113	fPIDmode(0),\r
114	\r
115	fNofTracks(0),\r
116	fPoolContent(0),\r
117	\r
118	fPhiMin(0),\r
119	fPhiMax(0),\r
120	\r
121	fPtTrigger(0),\r
122	fPhiTrigger(0),\r
123	fEtaTrigger(0),\r
124	\r
125	\r
126	fDeltaPhi(0),\r
127	fDeltaEta(0),\r
128	fk0InvMass(0)\r
129	{\r
130	fhadcuts = cuts; \r
131	}\r
132	\r
133	//_____________________________________________________\r
134	AliHFCorrelator::~AliHFCorrelator() \r
135	{\r
136	//\r
137	// destructor\r
138	// \r
139	\r
140	if(fPoolMgr) {delete fPoolMgr; fPoolMgr=0;} \r
141	if(fPool) {delete fPool; fPool=0;}\r
142	if(fhadcuts) {delete fhadcuts; fhadcuts=0;}\r
143	if(fAODEvent) {delete fAODEvent; fAODEvent=0;}\r
144	if(fAssociatedTracks) {delete fAssociatedTracks; fAssociatedTracks=0;}\r
145	if(fmcArray) {delete fmcArray; fmcArray=0;}\r
146	if(fReducedPart) {delete fReducedPart; fReducedPart=0;}\r
147	if(fD0cand) {delete fD0cand; fD0cand=0;}\r
148	\r
149	\r
150	if(fNofTracks) fNofTracks = 0;\r
151	\r
152	if(fPhiMin) fPhiMin = 0;\r
153	if(fPhiMax) fPhiMax = 0;\r
154	\r
155	if(fPtTrigger) fPtTrigger=0;\r
156	if(fPhiTrigger) fPhiTrigger=0;\r
157	if(fEtaTrigger) fEtaTrigger=0;\r
158	\r
159	if(fDeltaPhi) fDeltaPhi=0;\r
160	if(fDeltaEta) fDeltaEta=0;\r
161	\r
162	if(fk0InvMass) fk0InvMass=0;\r
163	\r
164	}\r
165	//_____________________________________________________\r
166	Bool_t AliHFCorrelator::DefineEventPool(){\r
167	// definition of the Pool Manager for Event Mixing\r
168	\r
169	\r
170	Int_t MaxNofEvents = fhadcuts->GetMaxNEventsInPool();\r
171	Int_t MinNofTracks = fhadcuts->GetMinNTracksInPool();\r
172	Int_t NofCentBins = fhadcuts->GetNCentPoolBins();\r
173	Double_t * CentBins = fhadcuts->GetCentPoolBins();\r
174	Int_t NofZVrtxBins = fhadcuts->GetNZvtxPoolBins();\r
175	Double_t *ZVrtxBins = fhadcuts->GetZvtxPoolBins();\r
176	\r
177	\r
178	fPoolMgr = new AliEventPoolManager(MaxNofEvents, MinNofTracks, NofCentBins, CentBins, NofZVrtxBins, ZVrtxBins);\r
179	if(!fPoolMgr) return kFALSE;\r
180	return kTRUE;\r
181	}\r
182	//_____________________________________________________\r
183	Bool_t AliHFCorrelator::Initialize(){\r
184	\r
185	// std::cout << "AliHFCorrelator::Initialize"<< std::endl;\r
186	// AliInfo("AliHFCorrelator::Initialize") ;\r
187	if(!fAODEvent){\r
188	AliInfo("No AOD event") ;\r
189	return kFALSE;\r
190	}\r
191	//std::cout << "No AOD event" << std::endl;\r
192	\r
193	AliCentrality *centralityObj = 0;\r
194	//Int_t multiplicity = -1;\r
195	Double_t MultipOrCent = -1;\r
196	\r
197	// initialize the pool for event mixing\r
198	if(!fsystem){ // pp\r
199	//multiplicity = fAODEvent->GetNTracks();\r
200	MultipOrCent = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(fAODEvent,-1.,1.);\r
201	// MultipOrCent = multiplicity; // convert from Int_t to Double_t\r
202	// AliInfo(Form("Multiplicity is %f", MultipOrCent));\r
203	}\r
204	if(fsystem){ // PbPb\r
205	\r
206	centralityObj = fAODEvent->GetHeader()->GetCentralityP();\r
207	MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");\r
208	// AliInfo(Form("Centrality is %f", MultipOrCent));\r
209	}\r
210	\r
211	AliAODVertex *vtx = fAODEvent->GetPrimaryVertex();\r
212	Double_t zvertex = vtx->GetZ(); // zvertex\r
213	Double_t * CentBins = fhadcuts->GetCentPoolBins();\r
214	Double_t poolmin=CentBins[0];\r
215	Double_t poolmax=CentBins[fhadcuts->GetNCentPoolBins()];\r
216	\r
217	\r
218	if(TMath::Abs(zvertex)>=10 \|\| MultipOrCent>poolmax \|\| MultipOrCent < poolmin) {\r
219	if(!fsystem)AliInfo(Form("pp or pA Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,MultipOrCent));\r
220	if(fsystem) AliInfo(Form("PbPb Event with Zvertex = %.2f cm and centrality = %.1f out of pool bounds, SKIPPING",zvertex,MultipOrCent));\r
221	\r
222	return kFALSE;\r
223	}\r
224	\r
225	fPool = fPoolMgr->GetEventPool(MultipOrCent, zvertex);\r
226	\r
227	if (!fPool){\r
228	AliInfo(Form("No pool found for multiplicity = %f, zVtx = %f cm", MultipOrCent, zvertex));\r
229	return kFALSE;\r
230	}\r
231	//fPool->PrintInfo();\r
232	return kTRUE;\r
233	}\r
234	\r
235	//_____________________________________________________\r
236	Bool_t AliHFCorrelator::ProcessEventPool(){\r
237	// analysis on Mixed Events\r
238	//cout << "AliHFCorrelator::ProcessEventPool"<< endl;\r
239	if(!fmixing) return kFALSE;\r
240	if(!fPool->IsReady()) return kFALSE;\r
241	if(fPool->GetCurrentNEvents()<fhadcuts->GetMinEventsToMix()) return kFALSE;\r
242	// fPool->PrintInfo();\r
243	fPoolContent = fPool->GetCurrentNEvents();\r
244	\r
245	return kTRUE;\r
246	\r
247	}\r
248	\r
249	//_____________________________________________________\r
250	Bool_t AliHFCorrelator::ProcessAssociatedTracks(Int_t EventLoopIndex, const TObjArray* associatedTracks){\r
d6e62ab6	251	// associatedTracks is not deleted, it should be (if needed) deleted in the user task\r
	252	\r
	253	if(!fmixing){ // analysis on Single Event\r
	254	if(fAssociatedTracks){\r
	255	fAssociatedTracks->Delete();\r
	256	delete fAssociatedTracks;\r
	257	} \r
	258	if(fselect==kHadron \|\| fselect ==kKaon){\r
	259	fAssociatedTracks = AcceptAndReduceTracks(fAODEvent);\r
	260	fAssociatedTracks->SetOwner(kTRUE);\r
	261	}\r
	262	if(fselect==kKZero) {\r
	263	fAssociatedTracks = AcceptAndReduceKZero(fAODEvent);\r
	264	fAssociatedTracks->SetOwner(kTRUE);\r
	265	} \r
	266	if(fselect==kElectron && associatedTracks) {\r
	267	fAssociatedTracks=new TObjArray(*associatedTracks);// Maybe better to call the copy constructor\r
	268	fAssociatedTracks->SetOwner(kFALSE);\r
	269	}\r
	270	\r
	271	}\r
	272	\r
	273	if(fmixing) { // analysis on Mixed Events\r
e9b1e63f	274	\r
e9b1e63f	275	\r
d6e62ab6	276	fAssociatedTracks = fPool->GetEvent(EventLoopIndex);\r
e9b1e63f	277	\r
d6e62ab6	278	\r
	279	\r
	280	\r
	281	} // end if mixing\r
	282	\r
	283	if(!fAssociatedTracks) return kFALSE;\r
	284	\r
	285	fNofTracks = fAssociatedTracks->GetEntriesFast(); \r
	286	\r
	287	return kTRUE;\r
e9b1e63f	288	\r
	289	}\r
	290	//_____________________________________________________\r
	291	Bool_t AliHFCorrelator::Correlate(Int_t loopindex){\r
	292	\r
	293	if(loopindex >= fNofTracks) return kFALSE;\r
	294	if(!fAssociatedTracks) return kFALSE;\r
	295	\r
	296	fReducedPart = (AliReducedParticle*)fAssociatedTracks->At(loopindex);\r
	297	\r
	298	\r
	299	fDeltaPhi = SetCorrectPhiRange(fPhiTrigger - fReducedPart->Phi());\r
	300	\r
	301	fDeltaEta = fEtaTrigger - fReducedPart->Eta();\r
	302	\r
	303	return kTRUE;\r
	304	\r
	305	}\r
	306	\r
	307	//_____________________________________________________\r
	308	Bool_t AliHFCorrelator::PoolUpdate(const TObjArray* associatedTracks){\r
	309	\r
	310	if(!fmixing) return kFALSE;\r
	311	if(!fPool) return kFALSE;\r
	312	if(fmixing) { // update the pool for Event Mixing\r
	313	TObjArray* objArr = NULL;\r
	314	if(fselect==kHadron \|\| fselect==kKaon) objArr = (TObjArray*)AcceptAndReduceTracks(fAODEvent);\r
	315	else if(fselect==kKZero) objArr = (TObjArray*)AcceptAndReduceKZero(fAODEvent);\r
	316	else if(fselect==kElectron && associatedTracks){\r
	317	objArr = new TObjArray(*associatedTracks);\r
	318	}\r
	319	else return kFALSE;\r
	320	if(objArr->GetEntriesFast()>0) fPool->UpdatePool(objArr); // updating the pool only if there are entries in the array\r
	321	}\r
	322	\r
	323	return kTRUE;\r
	324	\r
	325	}\r
	326	\r
	327	//_____________________________________________________\r
	328	Double_t AliHFCorrelator::SetCorrectPhiRange(Double_t phi){\r
	329	Double_t pi = TMath::Pi();\r
	330	\r
	331	if(phi<fPhiMin) phi = phi + 2*pi;\r
	332	if(phi>fPhiMax) phi = phi - 2*pi;\r
	333	\r
	334	return phi;\r
	335	}\r
	336	\r
	337	//_____________________________________________________\r
	338	TObjArray* AliHFCorrelator::AcceptAndReduceTracks(AliAODEvent* inputEvent){\r
	339	\r
	340	Double_t weight=1.;\r
	341	Int_t nTracks = inputEvent->GetNTracks();\r
	342	AliAODVertex * vtx = inputEvent->GetPrimaryVertex();\r
	343	Double_t pos[3],cov[6];\r
	344	vtx->GetXYZ(pos);\r
	345	vtx->GetCovarianceMatrix(cov);\r
	346	const AliESDVertex vESD(pos,cov,100.,100);\r
	347	\r
	348	Double_t Bz = inputEvent->GetMagneticField();\r
	349	\r
	350	\r
	351	TObjArray* tracksClone = new TObjArray;\r
352	tracksClone->SetOwner(kTRUE);\r
353	\r
354	//*******************************************************\r
355	// use reconstruction\r
356	if(fUseReco){\r
357	for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {\r
358	AliAODTrack* track = inputEvent->GetTrack(iTrack);\r
359	if (!track) continue;\r
360	if(!fhadcuts->IsHadronSelected(track,&vESD,Bz)) continue; // apply ESD level selections\r
361	if(!fhadcuts->Charge(fDCharge,track)) continue; // apply selection on charge, if required\r
362	\r
363	Double_t pT = track->Pt();\r
364	\r
365	//compute impact parameter\r
366	Double_t d0z0[2],covd0z0[3];\r
367	Double_t d0=-999999.;\r
368	if(fUseImpactParameter) track->PropagateToDCA(vtx,Bz,100,d0z0,covd0z0);\r
369	else d0z0[0] = 1. ; // random number - be careful with the cuts you applied\r
370	\r
371	if(fUseImpactParameter==1) d0 = TMath::Abs(d0z0[0]); // use impact parameter\r
372	if(fUseImpactParameter==2) { // use impact parameter over resolution\r
373	if(TMath::Abs(covd0z0[0])>0.00000001) d0 = TMath::Abs(d0z0[0])/TMath::Sqrt(covd0z0[0]); \r
374	else d0 = -1.; // if the resoultion is Zero, rejects the track - to be on the safe side\r
375	\r
376	}\r
377	\r
378	if(fmontecarlo) {// THIS TO BE CHECKED\r
379	Int_t hadLabel = track->GetLabel();\r
380	if(hadLabel < 0) continue; \r
381	}\r
382	\r
383	if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts\r
384	Bool_t rejectsoftpi = kTRUE;// TO BE CHECKED: DO WE WANT IT TO kTRUE AS A DEFAULT?\r
385	if(fD0cand && !fmixing) rejectsoftpi = fhadcuts->InvMassDstarRejection(fD0cand,track,fhypD0); // TO BE CHECKED: WHY NOT FOR EM?\r
386	\r
387	\r
388	if(fselect ==kKaon){ \r
389	if(!fhadcuts->CheckKaonCompatibility(track,fmontecarlo,fmcArray,fPIDmode)) continue; // check if it is a Kaon - data and MC\r
390	}\r
391	weight=fhadcuts->GetTrackWeight(pT,track->Eta(),pos[2]);\r
392	tracksClone->Add(new AliReducedParticle(track->Eta(), track->Phi(), pT,track->GetLabel(),track->GetID(),d0,rejectsoftpi,track->Charge(),weight));\r
393	} // end loop on tracks\r
394	} // end if use reconstruction kTRUE\r
395	\r
396	//*******************************************************\r
397	\r
398	//use MC truth\r
399	if(fmontecarlo && !fUseReco){\r
400	\r
401	for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) { \r
402	AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));\r
403	if (!mcPart) {\r
404	AliWarning("MC Particle not found in tree, skipping"); \r
405	continue;\r
406	}\r
407	if(!mcPart->Charge()) continue; // consider only charged tracks\r
408	\r
409	Int_t PDG =TMath::Abs(mcPart->PdgCode()); \r
410	if(fselect ==kHadron) {if(!((PDG==321)\|\|(PDG==211)\|\|(PDG==2212)\|\|(PDG==13)\|\|(PDG==11))) continue;} // select only if kaon, pion, proton, muon or electron\r
411	else if(fselect ==kKaon) {if(!(PDG==321)) continue;} // select only if kaon\r
412	else if(fselect ==kElectron) {if(!(PDG==11)) continue;} // select only if electron\r
413	\r
414	Double_t pT = mcPart->Pt();\r
415	Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet\r
416	if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts\r
417	\r
418	tracksClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,iPart,-1,d0,kFALSE,mcPart->Charge()));\r
419	}\r
420	\r
421	} // end if use MC truth\r
422	\r
423	\r
424	return tracksClone;\r
425	}\r
426	\r
427	//_____________________________________________________\r
428	TObjArray* AliHFCorrelator::AcceptAndReduceKZero(AliAODEvent* inputEvent){\r
429	\r
430	Int_t nOfVZeros = inputEvent->GetNumberOfV0s();\r
431	TObjArray* KZeroClone = new TObjArray;\r
432	AliAODVertex vertex1 = (AliAODVertex)inputEvent->GetPrimaryVertex();\r
433	\r
434	// use reconstruction \r
435	if(fUseReco){\r
436	Int_t v0label = -1;\r
437	Int_t pdgDgK0toPipi[2] = {211,211};\r
438	Double_t mPDGK0=0.497614;//TDatabasePDG::Instance()->GetParticle(310)->Mass();\r
439	const Int_t size = inputEvent->GetNumberOfV0s();\r
440	Int_t prevnegID[size];\r
441	Int_t prevposID[size];\r
442	for(Int_t iv0 =0; iv0< nOfVZeros; iv0++){// loop on all v0 candidates\r
443	AliAODv0 v0 = (static_cast<AliAODEvent> (inputEvent))->GetV0(iv0);\r
444	if(!v0) {\r
445	AliInfo(Form("is not a v0 at step, %d", iv0)) ;\r
446	//cout << "is not a v0 at step " << iv0 << endl;\r
447	continue;\r
448	}\r
449	\r
450	if(!fhadcuts->IsKZeroSelected(v0,vertex1)) continue; // check if it is a k0\r
451	\r
452	// checks to avoid double counting\r
453	Int_t negID = -999;\r
454	Int_t posID = -998;\r
455	//Int_t a = 0;\r
456	prevnegID[iv0] = -997;\r
457	prevposID[iv0]= -996;\r
458	negID = v0->GetNegID();\r
459	posID = v0->GetPosID();\r
460	Bool_t DoubleCounting = kFALSE;\r
461	\r
462	for(Int_t k=0; k<iv0; k++){\r
463	if(((negID==prevnegID[k])&&(posID==prevposID[k]))\|\|((negID==prevposID[k])&&(posID==prevnegID[k]))){\r
464	DoubleCounting = kTRUE;\r
465	//a=k;\r
466	break;\r
467	}//end if\r
468	} // end for\r
469	\r
470	if(DoubleCounting) continue;\r
471	else{ \r
472	prevposID[iv0] = posID; \r
473	prevnegID[iv0] = negID;\r
474	}\r
475	\r
476	if(fmontecarlo) v0label = v0->MatchToMC(310,fmcArray, 0, pdgDgK0toPipi); //match a K0 short\r
477	Double_t k0pt = v0->Pt();\r
478	Double_t k0eta = v0->Eta();\r
479	Double_t k0Phi = v0->Phi();\r
480	fk0InvMass = v0->MassK0Short(); \r
481	\r
482	//if (loopindex == 0) {\r
483	// if(!plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectra"))->Fill(k0InvMass,k0pt); // spectra for all k0\r
484	// if(plotassociation) ((TH2F)fOutput->FindObject("KZeroSpectraifHF"))->Fill(k0InvMass,k0pt); // spectra for k0 in association with a D\r
485	//}\r
486	// if there are more D* candidates per event, loopindex == 0 makes sure you fill the mass spectra only once!\r
487	\r
488	if(TMath::Abs(fk0InvMass-mPDGK0)>3*0.004) continue; // select candidates within 3 sigma\r
489	KZeroClone->Add(new AliReducedParticle(k0eta,k0Phi,k0pt,v0label));\r
490	\r
491	}\r
492	} // end if use reconstruction kTRUE\r
493	\r
494	\r
495	\r
496	//*********************************************************************//\r
497	//use MC truth\r
498	if(fmontecarlo && !fUseReco){\r
499	\r
500	for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) { \r
501	AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));\r
502	if (!mcPart) {\r
503	AliWarning("MC Particle not found in tree, skipping"); \r
504	continue;\r
505	}\r
506	\r
507	Int_t PDG =TMath::Abs(mcPart->PdgCode()); \r
508	if(!(PDG==310)) continue; // select only if k0 short\r
509	\r
510	Double_t pT = mcPart->Pt();\r
511	Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet\r
512	if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts\r
513	\r
514	KZeroClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,iPart,-1,d0,kTRUE,mcPart->Charge()));\r
515	}\r
516	\r
517	} // end if use MC truth\r
518	\r
519	\r
520	\r
521	return KZeroClone;\r
522	}\r