[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 = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(fAODEvent,-1.,1.);
	//	MultipOrCent = multiplicity; // convert from Int_t to Double_t
     //   AliInfo(Form("Multiplicity is %f", MultipOrCent));
	}
	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 or pA 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){
  // associatedTracks is not deleted, it should be (if needed) deleted in the user task
  
  if(!fmixing){ // analysis on Single Event
    if(fAssociatedTracks){
      fAssociatedTracks->Delete();
      delete fAssociatedTracks;
    }      
    if(fselect==kHadron || fselect ==kKaon){
      fAssociatedTracks = AcceptAndReduceTracks(fAODEvent);
      fAssociatedTracks->SetOwner(kTRUE);
    }
    if(fselect==kKZero) {
      fAssociatedTracks = AcceptAndReduceKZero(fAODEvent);
      fAssociatedTracks->SetOwner(kTRUE);
    }	
    if(fselect==kElectron && associatedTracks) {
      fAssociatedTracks=new TObjArray(*associatedTracks);// Maybe better to call the copy constructor
      fAssociatedTracks->SetOwner(kFALSE);
    }
    
  }
  
  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
896d3200	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"
	41	#include "AliAODMCParticle.h"
	42
	43	using std::cout;
	44	using std::endl;
	45
	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),
	61	fDCharge(0),
	62
	63	fmixing(kFALSE),
	64	fmontecarlo(kFALSE),
65	fsystem(kFALSE),
66	fUseReco(kTRUE),
67	fselect(kUndefined),
68
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),
105	fDCharge(0),
106
107	fmixing(kFALSE),
108	fmontecarlo(kFALSE),
109	fsystem(ppOrPbPb),
110	fUseReco(kTRUE),
111	fselect(kUndefined),
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
185	// std::cout << "AliHFCorrelator::Initialize"<< std::endl;
186	// AliInfo("AliHFCorrelator::Initialize") ;
187	if(!fAODEvent){
188	AliInfo("No AOD event") ;
189	return kFALSE;
190	}
191	//std::cout << "No AOD event" << std::endl;
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 = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(fAODEvent,-1.,1.);
201	// MultipOrCent = multiplicity; // convert from Int_t to Double_t
202	// AliInfo(Form("Multiplicity is %f", MultipOrCent));
203	}
204	if(fsystem){ // PbPb
205
206	centralityObj = fAODEvent->GetHeader()->GetCentralityP();
207	MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
208	// AliInfo(Form("Centrality is %f", MultipOrCent));
209	}
210
211	AliAODVertex *vtx = fAODEvent->GetPrimaryVertex();
212	Double_t zvertex = vtx->GetZ(); // zvertex
213	Double_t * CentBins = fhadcuts->GetCentPoolBins();
214	Double_t poolmin=CentBins[0];
215	Double_t poolmax=CentBins[fhadcuts->GetNCentPoolBins()];
216
217
218	if(TMath::Abs(zvertex)>=10 \|\| MultipOrCent>poolmax \|\| MultipOrCent < poolmin) {
219	if(!fsystem)AliInfo(Form("pp or pA Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,MultipOrCent));
220	if(fsystem) AliInfo(Form("PbPb Event with Zvertex = %.2f cm and centrality = %.1f out of pool bounds, SKIPPING",zvertex,MultipOrCent));
221
222	return kFALSE;
223	}
224
225	fPool = fPoolMgr->GetEventPool(MultipOrCent, zvertex);
226
227	if (!fPool){
228	AliInfo(Form("No pool found for multiplicity = %f, zVtx = %f cm", MultipOrCent, zvertex));
229	return kFALSE;
230	}
231	//fPool->PrintInfo();
232	return kTRUE;
233	}
234
235	//_____________________________________________________
236	Bool_t AliHFCorrelator::ProcessEventPool(){
237	// analysis on Mixed Events
238	//cout << "AliHFCorrelator::ProcessEventPool"<< endl;
239	if(!fmixing) return kFALSE;
240	if(!fPool->IsReady()) return kFALSE;
241	if(fPool->GetCurrentNEvents()<fhadcuts->GetMinEventsToMix()) return kFALSE;
242	// fPool->PrintInfo();
243	fPoolContent = fPool->GetCurrentNEvents();
244
245	return kTRUE;
246
247	}
248
249	//_____________________________________________________
250	Bool_t AliHFCorrelator::ProcessAssociatedTracks(Int_t EventLoopIndex, const TObjArray* associatedTracks){
251	// associatedTracks is not deleted, it should be (if needed) deleted in the user task
252
253	if(!fmixing){ // analysis on Single Event
254	if(fAssociatedTracks){
255	fAssociatedTracks->Delete();
256	delete fAssociatedTracks;
257	}
258	if(fselect==kHadron \|\| fselect ==kKaon){
259	fAssociatedTracks = AcceptAndReduceTracks(fAODEvent);
260	fAssociatedTracks->SetOwner(kTRUE);
261	}
262	if(fselect==kKZero) {
263	fAssociatedTracks = AcceptAndReduceKZero(fAODEvent);
264	fAssociatedTracks->SetOwner(kTRUE);
265	}
266	if(fselect==kElectron && associatedTracks) {
267	fAssociatedTracks=new TObjArray(*associatedTracks);// Maybe better to call the copy constructor
268	fAssociatedTracks->SetOwner(kFALSE);
269	}
270
271	}
272
273	if(fmixing) { // analysis on Mixed Events
274
275
276	fAssociatedTracks = fPool->GetEvent(EventLoopIndex);
277
278
279
280
281	} // end if mixing
282
283	if(!fAssociatedTracks) return kFALSE;
284
285	fNofTracks = fAssociatedTracks->GetEntriesFast();
286
287	return kTRUE;
288
289	}
290	//_____________________________________________________
291	Bool_t AliHFCorrelator::Correlate(Int_t loopindex){
292
293	if(loopindex >= fNofTracks) return kFALSE;
294	if(!fAssociatedTracks) return kFALSE;
295
296	fReducedPart = (AliReducedParticle*)fAssociatedTracks->At(loopindex);
297
298
299	fDeltaPhi = SetCorrectPhiRange(fPhiTrigger - fReducedPart->Phi());
300
301	fDeltaEta = fEtaTrigger - fReducedPart->Eta();
302
303	return kTRUE;
304
305	}
306
307	//_____________________________________________________
308	Bool_t AliHFCorrelator::PoolUpdate(const TObjArray* associatedTracks){
309
310	if(!fmixing) return kFALSE;
311	if(!fPool) return kFALSE;
312	if(fmixing) { // update the pool for Event Mixing
313	TObjArray* objArr = NULL;
314	if(fselect==kHadron \|\| fselect==kKaon) objArr = (TObjArray*)AcceptAndReduceTracks(fAODEvent);
315	else if(fselect==kKZero) objArr = (TObjArray*)AcceptAndReduceKZero(fAODEvent);
316	else if(fselect==kElectron && associatedTracks){
317	objArr = new TObjArray(*associatedTracks);
318	}
319	else return kFALSE;
320	if(objArr->GetEntriesFast()>0) fPool->UpdatePool(objArr); // updating the pool only if there are entries in the array
321	}
322
323	return kTRUE;
324
325	}
326
327	//_____________________________________________________
328	Double_t AliHFCorrelator::SetCorrectPhiRange(Double_t phi){
329	Double_t pi = TMath::Pi();
330
331	if(phi<fPhiMin) phi = phi + 2*pi;
332	if(phi>fPhiMax) phi = phi - 2*pi;
333
334	return phi;
335	}
336
337	//_____________________________________________________
338	TObjArray* AliHFCorrelator::AcceptAndReduceTracks(AliAODEvent* inputEvent){
339
340	Double_t weight=1.;
341	Int_t nTracks = inputEvent->GetNTracks();
342	AliAODVertex * vtx = inputEvent->GetPrimaryVertex();
343	Double_t pos[3],cov[6];
344	vtx->GetXYZ(pos);
345	vtx->GetCovarianceMatrix(cov);
346	const AliESDVertex vESD(pos,cov,100.,100);
347
348	Double_t Bz = inputEvent->GetMagneticField();
349
350
351	TObjArray* tracksClone = new TObjArray;
352	tracksClone->SetOwner(kTRUE);
353
354	//*******************************************************
355	// use reconstruction
356	if(fUseReco){
357	for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
358	AliAODTrack* track = inputEvent->GetTrack(iTrack);
359	if (!track) continue;
360	if(!fhadcuts->IsHadronSelected(track,&vESD,Bz)) continue; // apply ESD level selections
361	if(!fhadcuts->Charge(fDCharge,track)) continue; // apply selection on charge, if required
362
363	Double_t pT = track->Pt();
364
365	//compute impact parameter
366	Double_t d0z0[2],covd0z0[3];
367	Double_t d0=-999999.;
368	if(fUseImpactParameter) track->PropagateToDCA(vtx,Bz,100,d0z0,covd0z0);
369	else d0z0[0] = 1. ; // random number - be careful with the cuts you applied
370
371	if(fUseImpactParameter==1) d0 = TMath::Abs(d0z0[0]); // use impact parameter
372	if(fUseImpactParameter==2) { // use impact parameter over resolution
373	if(TMath::Abs(covd0z0[0])>0.00000001) d0 = TMath::Abs(d0z0[0])/TMath::Sqrt(covd0z0[0]);
374	else d0 = -1.; // if the resoultion is Zero, rejects the track - to be on the safe side
375
376	}
377
378	if(fmontecarlo) {// THIS TO BE CHECKED
379	Int_t hadLabel = track->GetLabel();
380	if(hadLabel < 0) continue;
381	}
382
383	if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
384	Bool_t rejectsoftpi = kTRUE;// TO BE CHECKED: DO WE WANT IT TO kTRUE AS A DEFAULT?
385	if(fD0cand && !fmixing) rejectsoftpi = fhadcuts->InvMassDstarRejection(fD0cand,track,fhypD0); // TO BE CHECKED: WHY NOT FOR EM?
386
387
388	if(fselect ==kKaon){
389	if(!fhadcuts->CheckKaonCompatibility(track,fmontecarlo,fmcArray,fPIDmode)) continue; // check if it is a Kaon - data and MC
390	}
391	weight=fhadcuts->GetTrackWeight(pT,track->Eta(),pos[2]);
392	tracksClone->Add(new AliReducedParticle(track->Eta(), track->Phi(), pT,track->GetLabel(),track->GetID(),d0,rejectsoftpi,track->Charge(),weight));
393	} // end loop on tracks
394	} // end if use reconstruction kTRUE
395
396	//*******************************************************
397
398	//use MC truth
399	if(fmontecarlo && !fUseReco){
400
401	for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) {
402	AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));
403	if (!mcPart) {
404	AliWarning("MC Particle not found in tree, skipping");
405	continue;
406	}
407	if(!mcPart->Charge()) continue; // consider only charged tracks
408
409	Int_t PDG =TMath::Abs(mcPart->PdgCode());
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
411	else if(fselect ==kKaon) {if(!(PDG==321)) continue;} // select only if kaon
412	else if(fselect ==kElectron) {if(!(PDG==11)) continue;} // select only if electron
413
414	Double_t pT = mcPart->Pt();
415	Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet
416	if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
417
418	tracksClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,iPart,-1,d0,kFALSE,mcPart->Charge()));
419	}
420
421	} // end if use MC truth
422
423
424	return tracksClone;
425	}
426
427	//_____________________________________________________
428	TObjArray* AliHFCorrelator::AcceptAndReduceKZero(AliAODEvent* inputEvent){
429
430	Int_t nOfVZeros = inputEvent->GetNumberOfV0s();
431	TObjArray* KZeroClone = new TObjArray;
432	AliAODVertex vertex1 = (AliAODVertex)inputEvent->GetPrimaryVertex();
433
434	// use reconstruction
435	if(fUseReco){
436	Int_t v0label = -1;
437	Int_t pdgDgK0toPipi[2] = {211,211};
438	Double_t mPDGK0=0.497614;//TDatabasePDG::Instance()->GetParticle(310)->Mass();
439	const Int_t size = inputEvent->GetNumberOfV0s();
440	Int_t prevnegID[size];
441	Int_t prevposID[size];
442	for(Int_t iv0 =0; iv0< nOfVZeros; iv0++){// loop on all v0 candidates
443	AliAODv0 v0 = (static_cast<AliAODEvent> (inputEvent))->GetV0(iv0);
444	if(!v0) {
445	AliInfo(Form("is not a v0 at step, %d", iv0)) ;
446	//cout << "is not a v0 at step " << iv0 << endl;
447	continue;
448	}
449
450	if(!fhadcuts->IsKZeroSelected(v0,vertex1)) continue; // check if it is a k0
451
452	// checks to avoid double counting
453	Int_t negID = -999;
454	Int_t posID = -998;
455	//Int_t a = 0;
456	prevnegID[iv0] = -997;
457	prevposID[iv0]= -996;
458	negID = v0->GetNegID();
459	posID = v0->GetPosID();
460	Bool_t DoubleCounting = kFALSE;
461
462	for(Int_t k=0; k<iv0; k++){
463	if(((negID==prevnegID[k])&&(posID==prevposID[k]))\|\|((negID==prevposID[k])&&(posID==prevnegID[k]))){
464	DoubleCounting = kTRUE;
465	//a=k;
466	break;
467	}//end if
468	} // end for
469
470	if(DoubleCounting) continue;
471	else{
472	prevposID[iv0] = posID;
473	prevnegID[iv0] = negID;
474	}
475
476	if(fmontecarlo) v0label = v0->MatchToMC(310,fmcArray, 0, pdgDgK0toPipi); //match a K0 short
477	Double_t k0pt = v0->Pt();
478	Double_t k0eta = v0->Eta();
479	Double_t k0Phi = v0->Phi();
480	fk0InvMass = v0->MassK0Short();
481
482	//if (loopindex == 0) {
483	// if(!plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectra"))->Fill(k0InvMass,k0pt); // spectra for all k0
484	// if(plotassociation) ((TH2F)fOutput->FindObject("KZeroSpectraifHF"))->Fill(k0InvMass,k0pt); // spectra for k0 in association with a D
485	//}
486	// if there are more D* candidates per event, loopindex == 0 makes sure you fill the mass spectra only once!
487
488	if(TMath::Abs(fk0InvMass-mPDGK0)>3*0.004) continue; // select candidates within 3 sigma
489	KZeroClone->Add(new AliReducedParticle(k0eta,k0Phi,k0pt,v0label));
490
491	}
492	} // end if use reconstruction kTRUE
493
494
495
496	//*********************************************************************//
497	//use MC truth
498	if(fmontecarlo && !fUseReco){
499
500	for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) {
501	AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));
502	if (!mcPart) {
503	AliWarning("MC Particle not found in tree, skipping");
504	continue;
505	}
506
507	Int_t PDG =TMath::Abs(mcPart->PdgCode());
508	if(!(PDG==310)) continue; // select only if k0 short
509
510	Double_t pT = mcPart->Pt();
511	Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet
512	if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
513
514	KZeroClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,iPart,-1,d0,kTRUE,mcPart->Charge()));
515	}
516
517	} // end if use MC truth
518
519
520
521	return KZeroClone;
522	}