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

/**************************************************************************
 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
//
//
//             Base class for Heavy Flavour Correlations Analysis
//             Single Event and Mixed Event Analysis are implemented
//
//-----------------------------------------------------------------------
//          
//
//						   Author S.Bjelogrlic
//                         Utrecht University 
//                      sandro.bjelogrlic@cern.ch
//
//-----------------------------------------------------------------------

/* $Id$ */

#include <TParticle.h>
#include <TVector3.h>
#include <TChain.h>
#include "TROOT.h"
#include "AliHFCorrelator.h"
#include "AliRDHFCutsDStartoKpipi.h"
#include "AliHFAssociatedTrackCuts.h"
#include "AliEventPoolManager.h"
#include "AliReducedParticle.h"
#include "AliCentrality.h"
#include "AliAODMCParticle.h"

using std::cout;
using std::endl;

//_____________________________________________________
AliHFCorrelator::AliHFCorrelator() :
//
// default constructor
//
TNamed(),
fPoolMgr(0x0),         
fPool(0x0),
fhadcuts(0x0),
fAODEvent(0x0),
fAssociatedTracks(0x0),
fmcArray(0x0),
fReducedPart(0x0),
fD0cand(0x0), 
fhypD0(0), 
fDCharge(0),

fmixing(kFALSE),
fmontecarlo(kFALSE),
fsystem(kFALSE),
fUseReco(kTRUE),
fselect(kUndefined),

fUseImpactParameter(0),
fPIDmode(0),

fNofTracks(0),
fPoolContent(0),

fPhiMin(0),
fPhiMax(0),

fPtTrigger(0),
fPhiTrigger(0),
fEtaTrigger(0),


fDeltaPhi(0),
fDeltaEta(0),
fk0InvMass(0)

{
	// default constructor	
}


//_____________________________________________________
AliHFCorrelator::AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t ppOrPbPb) :
TNamed(name,"title"),
fPoolMgr(0x0),         
fPool(0x0),
fhadcuts(0x0),
fAODEvent(0x0),
fAssociatedTracks(0x0),
fmcArray(0x0),
fReducedPart(0x0),
fD0cand(0x0), 
fhypD0(0),
fDCharge(0),

fmixing(kFALSE),
fmontecarlo(kFALSE),
fsystem(ppOrPbPb),
fUseReco(kTRUE),
fselect(kUndefined),
fUseImpactParameter(0),
fPIDmode(0),

fNofTracks(0),
fPoolContent(0),

fPhiMin(0),
fPhiMax(0),

fPtTrigger(0),
fPhiTrigger(0),
fEtaTrigger(0),


fDeltaPhi(0),
fDeltaEta(0),
fk0InvMass(0)
{
	fhadcuts = cuts; 
}

//_____________________________________________________
AliHFCorrelator::~AliHFCorrelator() 
{
//
// destructor
//	
	
	if(fPoolMgr)  {delete fPoolMgr; fPoolMgr=0;}       
	if(fPool) {delete fPool; fPool=0;}
	if(fhadcuts) {delete fhadcuts; fhadcuts=0;}
	if(fAODEvent) {delete fAODEvent; fAODEvent=0;}
	if(fAssociatedTracks) {delete fAssociatedTracks; fAssociatedTracks=0;}
	if(fmcArray) {delete fmcArray; fmcArray=0;}
	if(fReducedPart) {delete fReducedPart; fReducedPart=0;}
	if(fD0cand) {delete fD0cand; fD0cand=0;}
	
	
	if(fNofTracks) fNofTracks = 0;
	
	if(fPhiMin) fPhiMin = 0;
	if(fPhiMax) fPhiMax = 0;
	
	if(fPtTrigger) fPtTrigger=0;
	if(fPhiTrigger) fPhiTrigger=0;
	if(fEtaTrigger) fEtaTrigger=0;
	
	if(fDeltaPhi) fDeltaPhi=0;
	if(fDeltaEta) fDeltaEta=0;
	
	if(fk0InvMass) fk0InvMass=0;
	
}
//_____________________________________________________
Bool_t AliHFCorrelator::DefineEventPool(){
	// definition of the Pool Manager for Event Mixing
	

	Int_t MaxNofEvents = fhadcuts->GetMaxNEventsInPool();
	Int_t MinNofTracks = fhadcuts->GetMinNTracksInPool();
	Int_t NofCentBins = fhadcuts->GetNCentPoolBins();
	Double_t * CentBins = fhadcuts->GetCentPoolBins();
	Int_t NofZVrtxBins = fhadcuts->GetNZvtxPoolBins();
	Double_t *ZVrtxBins = fhadcuts->GetZvtxPoolBins();
		
			
	fPoolMgr = new AliEventPoolManager(MaxNofEvents, MinNofTracks, NofCentBins, CentBins, NofZVrtxBins, ZVrtxBins);
	if(!fPoolMgr) return kFALSE;
	return kTRUE;
}
//_____________________________________________________
Bool_t AliHFCorrelator::Initialize(){
	
    //  std::cout << "AliHFCorrelator::Initialize"<< std::endl;
  AliInfo("AliHFCorrelator::Initialize") ;
  if(!fAODEvent){
    AliInfo("No AOD event") ;
    return kFALSE;
  }
    //std::cout << "No AOD event" << std::endl;
	
	AliCentrality *centralityObj = 0;
	Int_t multiplicity = -1;
	Double_t MultipOrCent = -1;
	
	// initialize the pool for event mixing
	if(!fsystem){ // pp
	multiplicity = fAODEvent->GetNTracks();
		MultipOrCent = multiplicity; // convert from Int_t to Double_t
	}
	if(fsystem){ // PbPb
		
		centralityObj = fAODEvent->GetHeader()->GetCentralityP();
		MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
		AliInfo(Form("Centrality is %f", MultipOrCent));
	}
	
	AliAODVertex *vtx = fAODEvent->GetPrimaryVertex();
	Double_t zvertex = vtx->GetZ(); // zvertex
	Double_t * CentBins = fhadcuts->GetCentPoolBins();
	Double_t poolmin=CentBins[0];
	Double_t poolmax=CentBins[fhadcuts->GetNCentPoolBins()];

	
		if(TMath::Abs(zvertex)>=10 || MultipOrCent>poolmax || MultipOrCent < poolmin) {
		if(!fsystem)AliInfo(Form("pp Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,MultipOrCent));
		if(fsystem) AliInfo(Form("PbPb Event with Zvertex = %.2f cm and centrality = %.1f  out of pool bounds, SKIPPING",zvertex,MultipOrCent));

			return kFALSE;
		}
	
	fPool = fPoolMgr->GetEventPool(MultipOrCent, zvertex);
	
	if (!fPool){
		AliInfo(Form("No pool found for multiplicity = %f, zVtx = %f cm", MultipOrCent, zvertex));
	    return kFALSE;
	}
	//fPool->PrintInfo();
	return kTRUE;
}

//_____________________________________________________
Bool_t AliHFCorrelator::ProcessEventPool(){
	 // analysis on Mixed Events
	//cout << "AliHFCorrelator::ProcessEventPool"<< endl;
		if(!fmixing) return kFALSE;
		if(!fPool->IsReady()) return kFALSE;
		if(fPool->GetCurrentNEvents()<fhadcuts->GetMinEventsToMix()) return kFALSE;
	//	fPool->PrintInfo();
		fPoolContent = fPool->GetCurrentNEvents();
		
		return kTRUE;
	
}

//_____________________________________________________
Bool_t AliHFCorrelator::ProcessAssociatedTracks(Int_t EventLoopIndex, const TObjArray* associatedTracks){
  // TODO: memory leak needs to be fixed, for every call, a new array
  // is allocated, but the pointer immediately lost. The cleanup is
  // not straightforward as in the case of event mixing the pointer
  // will be an external array which must not be deleted.
	fAssociatedTracks = new TObjArray();

	if(!fmixing){ // analysis on Single Event
		
		
		if(fselect==kHadron || fselect ==kKaon)	fAssociatedTracks = AcceptAndReduceTracks(fAODEvent);
		if(fselect==kKZero) {fAssociatedTracks = AcceptAndReduceKZero(fAODEvent);}	
		if(fselect==kElectron && associatedTracks) fAssociatedTracks=new TObjArray(*associatedTracks);
		
	}
	
	if(fmixing) { // analysis on Mixed Events
		
			
		fAssociatedTracks = fPool->GetEvent(EventLoopIndex);
				
				
	} // end if mixing
	
	if(!fAssociatedTracks) return kFALSE;
	
	fNofTracks = fAssociatedTracks->GetEntriesFast(); 
		
	return kTRUE;
	
}
//_____________________________________________________
Bool_t AliHFCorrelator::Correlate(Int_t loopindex){

	if(loopindex >= fNofTracks) return kFALSE;
	if(!fAssociatedTracks) return kFALSE;
	
	fReducedPart = (AliReducedParticle*)fAssociatedTracks->At(loopindex);
	

	fDeltaPhi = SetCorrectPhiRange(fPhiTrigger - fReducedPart->Phi());
	
	fDeltaEta = fEtaTrigger - fReducedPart->Eta();

	return kTRUE;
	
}
		
//_____________________________________________________
Bool_t AliHFCorrelator::PoolUpdate(const TObjArray* associatedTracks){

	if(!fmixing) return kFALSE;
	if(!fPool) return kFALSE;
	if(fmixing) { // update the pool for Event Mixing
		TObjArray* objArr = NULL;
		if(fselect==kHadron || fselect==kKaon) objArr = (TObjArray*)AcceptAndReduceTracks(fAODEvent);
		else if(fselect==kKZero) objArr = (TObjArray*)AcceptAndReduceKZero(fAODEvent);
		else if(fselect==kElectron && associatedTracks){
		  objArr = new TObjArray(*associatedTracks);
		}
		else return kFALSE;
		if(objArr->GetEntriesFast()>0) fPool->UpdatePool(objArr); // updating the pool only if there are entries in the array
	}
		
	return kTRUE;
	
}
		
//_____________________________________________________
Double_t AliHFCorrelator::SetCorrectPhiRange(Double_t phi){
	Double_t pi = TMath::Pi();
	
	if(phi<fPhiMin) phi = phi + 2*pi;
	if(phi>fPhiMax) phi = phi - 2*pi;
	
	return phi;
}

//_____________________________________________________
TObjArray*  AliHFCorrelator::AcceptAndReduceTracks(AliAODEvent* inputEvent){

  Double_t weight=1.;
  Int_t nTracks = inputEvent->GetNTracks();
  AliAODVertex * vtx = inputEvent->GetPrimaryVertex();
  Double_t pos[3],cov[6];
  vtx->GetXYZ(pos);
  vtx->GetCovarianceMatrix(cov);
  const AliESDVertex vESD(pos,cov,100.,100);
  
  Double_t Bz = inputEvent->GetMagneticField();
	
  
  TObjArray* tracksClone = new TObjArray;
  tracksClone->SetOwner(kTRUE);
  
  //*******************************************************
  // use reconstruction
  if(fUseReco){
    for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
      AliAODTrack* track = inputEvent->GetTrack(iTrack);
      if (!track) continue;
      if(!fhadcuts->IsHadronSelected(track,&vESD,Bz)) continue; // apply ESD level selections
      if(!fhadcuts->Charge(fDCharge,track)) continue; // apply selection on charge, if required

      Double_t pT = track->Pt();
      
      //compute impact parameter
      Double_t d0z0[2],covd0z0[3];
      Double_t d0=-999999.;
      if(fUseImpactParameter) track->PropagateToDCA(vtx,Bz,100,d0z0,covd0z0);
      else d0z0[0] = 1. ; // random number - be careful with the cuts you applied
      
      if(fUseImpactParameter==1) d0 = TMath::Abs(d0z0[0]); // use impact parameter
      if(fUseImpactParameter==2) { // use impact parameter over resolution
	if(TMath::Abs(covd0z0[0])>0.00000001) d0 = TMath::Abs(d0z0[0])/TMath::Sqrt(covd0z0[0]); 
	else d0 = -1.; // if the resoultion is Zero, rejects the track - to be on the safe side
	
      }
      
      if(fmontecarlo) {// THIS TO BE CHECKED
        Int_t hadLabel = track->GetLabel();
	if(hadLabel < 0) continue;	
      }
      
      if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
      Bool_t rejectsoftpi = kTRUE;// TO BE CHECKED: DO WE WANT IT TO kTRUE AS A DEFAULT?
      if(fD0cand && !fmixing) rejectsoftpi = fhadcuts->InvMassDstarRejection(fD0cand,track,fhypD0); // TO BE CHECKED: WHY NOT FOR EM?
      
      
      if(fselect ==kKaon){	
	if(!fhadcuts->CheckKaonCompatibility(track,fmontecarlo,fmcArray,fPIDmode)) continue; // check if it is a Kaon - data and MC
      }
      weight=fhadcuts->GetTrackWeight(pT,track->Eta(),pos[2]);
      tracksClone->Add(new AliReducedParticle(track->Eta(), track->Phi(), pT,track->GetLabel(),track->GetID(),d0,rejectsoftpi,track->Charge(),weight));
    } // end loop on tracks
  } // end if use reconstruction kTRUE
  
  //*******************************************************
  
  //use MC truth
  if(fmontecarlo && !fUseReco){
    
    for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) { 
      AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));
      if (!mcPart) {
	AliWarning("MC Particle not found in tree, skipping"); 
	continue;
      }
      if(!mcPart->Charge()) continue; // consider only charged tracks
      
      Int_t PDG =TMath::Abs(mcPart->PdgCode()); 
if(fselect ==kHadron) {if(!((PDG==321)||(PDG==211)||(PDG==2212)||(PDG==13)||(PDG==11))) continue;} // select only if kaon, pion, proton, muon or electron
      else if(fselect ==kKaon) {if(!(PDG==321)) continue;} // select only if kaon
      else if(fselect ==kElectron) {if(!(PDG==11)) continue;} // select only if electron

      Double_t pT = mcPart->Pt();
      Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet
      if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
      
      tracksClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,iPart,-1,d0,kFALSE,mcPart->Charge()));
    }
    
  } // end if use  MC truth
  
  
  return tracksClone;
}

//_____________________________________________________
TObjArray*  AliHFCorrelator::AcceptAndReduceKZero(AliAODEvent* inputEvent){
	
	Int_t nOfVZeros = inputEvent->GetNumberOfV0s();
	TObjArray* KZeroClone = new TObjArray;
	AliAODVertex *vertex1 = (AliAODVertex*)inputEvent->GetPrimaryVertex();

 // use reconstruction	 	
     if(fUseReco){
	Int_t v0label = -1;
	Int_t pdgDgK0toPipi[2] = {211,211};
	Double_t mPDGK0=0.497614;//TDatabasePDG::Instance()->GetParticle(310)->Mass();
	const Int_t size = inputEvent->GetNumberOfV0s();
	Int_t prevnegID[size];
	Int_t prevposID[size];
	for(Int_t iv0 =0; iv0< nOfVZeros; iv0++){// loop on all v0 candidates
		AliAODv0 *v0 = (static_cast<AliAODEvent*> (inputEvent))->GetV0(iv0);
		if(!v0) {
		  AliInfo(Form("is not a v0 at step, %d", iv0)) ;
		  //cout << "is not a v0 at step " << iv0 << endl;
		  continue;
		}
		
		if(!fhadcuts->IsKZeroSelected(v0,vertex1)) continue; // check if it is a k0
	    
		// checks to avoid double counting
		Int_t negID = -999;
		Int_t posID = -998;
		//Int_t a = 0;
		prevnegID[iv0] = -997;
		prevposID[iv0]= -996;
		negID = v0->GetNegID();
		posID = v0->GetPosID();
		Bool_t DoubleCounting = kFALSE;
		
		for(Int_t k=0; k<iv0; k++){
			if(((negID==prevnegID[k])&&(posID==prevposID[k]))||((negID==prevposID[k])&&(posID==prevnegID[k]))){
				DoubleCounting = kTRUE;
				//a=k;
				break;
			}//end if
		} // end for
		
		if(DoubleCounting) continue;
		else{ 
			prevposID[iv0] = posID; 
			prevnegID[iv0] = negID;
		}
		
		if(fmontecarlo)	v0label = v0->MatchToMC(310,fmcArray, 0, pdgDgK0toPipi); //match a K0 short
		Double_t k0pt = v0->Pt();
		Double_t k0eta = v0->Eta();
		Double_t k0Phi = v0->Phi();
	    fk0InvMass = v0->MassK0Short();	
		
		//if (loopindex == 0) {
		//	if(!plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectra"))->Fill(k0InvMass,k0pt); // spectra for all k0
		//	if(plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectraifHF"))->Fill(k0InvMass,k0pt); // spectra for k0 in association with a D*
		//}
		// if there are more D* candidates per event, loopindex == 0 makes sure you fill the mass spectra only once!
		
		if(TMath::Abs(fk0InvMass-mPDGK0)>3*0.004) continue; // select candidates within 3 sigma
		KZeroClone->Add(new AliReducedParticle(k0eta,k0Phi,k0pt,v0label));
		
	}
     } // end if use reconstruction kTRUE
	

//*********************************************************************//
     //use MC truth
     if(fmontecarlo && !fUseReco){
		
		for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) { 
			AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));
			if (!mcPart) {
				AliWarning("MC Particle not found in tree, skipping"); 
				continue;
			}
			
			Int_t PDG =TMath::Abs(mcPart->PdgCode()); 
			if(!(PDG==310)) continue; // select only if k0 short
		
			Double_t pT = mcPart->Pt();
            Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet
			if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
			
			KZeroClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,iPart,-1,d0,kTRUE,mcPart->Charge()));
		}
		
	} // end if use  MC truth


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