[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: AliHFCorrelator.cxx 64115 2013-09-05 12:34:55Z arossi $ */

#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),
fDMesonCutObject(0x0),
fAssociatedTracks(0x0),
fmcArray(0x0),
fReducedPart(0x0),
fD0cand(0x0), 
fhypD0(0), 
fDCharge(0),

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

fUseImpactParameter(0),
fPIDmode(0),

fNofTracks(0),
fPoolContent(0),

fPhiMin(0),
fPhiMax(0),

fMultCentr(-1),

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


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

{
	// default constructor	
}


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

fmixing(kFALSE),
fmontecarlo(kFALSE),
fUseCentrality(useCentrality),
fUseReco(kTRUE),
fselect(kUndefined),
fUseImpactParameter(0),
fPIDmode(0),

fNofTracks(0),
fPoolContent(0),

fPhiMin(0),
fPhiMax(0),

fMultCentr(-1),

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


fDeltaPhi(0),
fDeltaEta(0),
fk0InvMass(0)
{
	fhadcuts = cuts;
     if(!fDMesonCutObject) AliInfo("D meson cut object not loaded - if using centrality the estimator will be V0M!");
}

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

fmixing(kFALSE),
fmontecarlo(kFALSE),
fUseCentrality(useCentrality),
fUseReco(kTRUE),
fselect(kUndefined),
fUseImpactParameter(0),
fPIDmode(0),

fNofTracks(0),
fPoolContent(0),

fPhiMin(0),
fPhiMax(0),

fMultCentr(-1),

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


fDeltaPhi(0),
fDeltaEta(0),
fk0InvMass(0)
{
	fhadcuts = cuts;
    fDMesonCutObject = cutObject;
    
    if(!fDMesonCutObject) AliInfo("D meson cut object not implemented properly! Check your centrality estimators");
}


//_____________________________________________________
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(fDMesonCutObject) {delete fDMesonCutObject; fDMesonCutObject=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(!fUseCentrality){ // pp, pA
	//multiplicity = fAODEvent->GetNumberOfTracks();
        //MultipOrCent = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(fAODEvent,-1.,1.);
        fMultCentr = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(fAODEvent,-1.,1.);
	//	MultipOrCent = multiplicity; // convert from Int_t to Double_t
     //   AliInfo(Form("Multiplicity is %f", MultipOrCent));
	}
	if(fUseCentrality){ // PbPb
		if(!fDMesonCutObject){
           
                centralityObj = ((AliVAODHeader*)fAODEvent->GetHeader())->GetCentralityP();
		fMultCentr = centralityObj->GetCentralityPercentileUnchecked("V0M");
        }
        else fMultCentr = fDMesonCutObject->GetCentrality(fAODEvent);
//		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 || fMultCentr>poolmax || fMultCentr < poolmin) {
		if(!fUseCentrality)AliInfo(Form("Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,fMultCentr));
		if(fUseCentrality) AliInfo(Form("Event with Zvertex = %.2f cm and centrality = %.1f  out of pool bounds, SKIPPING",zvertex,fMultCentr));

			return kFALSE;
		}
	
	fPool = fPoolMgr->GetEventPool(fMultCentr, zvertex);
	
	if (!fPool){
		AliInfo(Form("No pool found for multiplicity = %f, zVtx = %f cm", fMultCentr, 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->GetNumberOfTracks();
  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 = dynamic_cast<AliAODTrack*>(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
c683985a	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
1b8b4a2c	29	/* $Id: AliHFCorrelator.cxx 64115 2013-09-05 12:34:55Z arossi $ */
c683985a	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),
1b8b4a2c	56	fDMesonCutObject(0x0),
c683985a	57	fAssociatedTracks(0x0),
	58	fmcArray(0x0),
	59	fReducedPart(0x0),
	60	fD0cand(0x0),
	61	fhypD0(0),
	62	fDCharge(0),
	63
	64	fmixing(kFALSE),
	65	fmontecarlo(kFALSE),
1b8b4a2c	66	fUseCentrality(kFALSE),
c683985a	67	fUseReco(kTRUE),
	68	fselect(kUndefined),
	69
	70	fUseImpactParameter(0),
	71	fPIDmode(0),
	72
	73	fNofTracks(0),
	74	fPoolContent(0),
	75
	76	fPhiMin(0),
	77	fPhiMax(0),
	78
1b8b4a2c EB	79	fMultCentr(-1),
1b8b4a2c EB	80
c683985a	81	fPtTrigger(0),
	82	fPhiTrigger(0),
	83	fEtaTrigger(0),
	84
	85
	86	fDeltaPhi(0),
	87	fDeltaEta(0),
	88	fk0InvMass(0)
	89
	90	{
	91	// default constructor
	92	}
	93
	94
	95
	96	//_____________________________________________________
1b8b4a2c	97	AliHFCorrelator::AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t useCentrality) :
c683985a	98	TNamed(name,"title"),
	99	fPoolMgr(0x0),
	100	fPool(0x0),
	101	fhadcuts(0x0),
	102	fAODEvent(0x0),
1b8b4a2c	103	fDMesonCutObject(0x0),
c683985a	104	fAssociatedTracks(0x0),
	105	fmcArray(0x0),
	106	fReducedPart(0x0),
	107	fD0cand(0x0),
	108	fhypD0(0),
	109	fDCharge(0),
	110
	111	fmixing(kFALSE),
	112	fmontecarlo(kFALSE),
1b8b4a2c EB	113	fUseCentrality(useCentrality),
	114	fUseReco(kTRUE),
	115	fselect(kUndefined),
	116	fUseImpactParameter(0),
	117	fPIDmode(0),
	118
	119	fNofTracks(0),
	120	fPoolContent(0),
	121
	122	fPhiMin(0),
	123	fPhiMax(0),
	124
	125	fMultCentr(-1),
	126
	127	fPtTrigger(0),
	128	fPhiTrigger(0),
	129	fEtaTrigger(0),
	130
	131
	132	fDeltaPhi(0),
	133	fDeltaEta(0),
	134	fk0InvMass(0)
	135	{
	136	fhadcuts = cuts;
	137	if(!fDMesonCutObject) AliInfo("D meson cut object not loaded - if using centrality the estimator will be V0M!");
	138	}
	139
	140	//_______________________________________________________________________________________
	141	AliHFCorrelator::AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts cuts, Bool_t useCentrality, AliRDHFCuts cutObject) :
	142	TNamed(name,"title"),
	143	fPoolMgr(0x0),
	144	fPool(0x0),
	145	fhadcuts(0x0),
	146	fAODEvent(0x0),
	147	fDMesonCutObject(0x0),
	148	fAssociatedTracks(0x0),
	149	fmcArray(0x0),
	150	fReducedPart(0x0),
	151	fD0cand(0x0),
	152	fhypD0(0),
	153	fDCharge(0),
	154
	155	fmixing(kFALSE),
	156	fmontecarlo(kFALSE),
	157	fUseCentrality(useCentrality),
c683985a	158	fUseReco(kTRUE),
	159	fselect(kUndefined),
	160	fUseImpactParameter(0),
	161	fPIDmode(0),
	162
	163	fNofTracks(0),
	164	fPoolContent(0),
	165
	166	fPhiMin(0),
	167	fPhiMax(0),
	168
1b8b4a2c EB	169	fMultCentr(-1),
1b8b4a2c EB	170
c683985a	171	fPtTrigger(0),
	172	fPhiTrigger(0),
	173	fEtaTrigger(0),
	174
	175
	176	fDeltaPhi(0),
	177	fDeltaEta(0),
	178	fk0InvMass(0)
	179	{
1b8b4a2c EB	180	fhadcuts = cuts;
	181	fDMesonCutObject = cutObject;
	182
	183	if(!fDMesonCutObject) AliInfo("D meson cut object not implemented properly! Check your centrality estimators");
c683985a	184	}
c683985a	185
1b8b4a2c EB	186
1b8b4a2c EB	187
c683985a	188	//_____________________________________________________
	189	AliHFCorrelator::~AliHFCorrelator()
	190	{
	191	//
	192	// destructor
	193	//
	194
	195	if(fPoolMgr) {delete fPoolMgr; fPoolMgr=0;}
	196	if(fPool) {delete fPool; fPool=0;}
	197	if(fhadcuts) {delete fhadcuts; fhadcuts=0;}
	198	if(fAODEvent) {delete fAODEvent; fAODEvent=0;}
1b8b4a2c	199	if(fDMesonCutObject) {delete fDMesonCutObject; fDMesonCutObject=0;}
c683985a	200	if(fAssociatedTracks) {delete fAssociatedTracks; fAssociatedTracks=0;}
	201	if(fmcArray) {delete fmcArray; fmcArray=0;}
	202	if(fReducedPart) {delete fReducedPart; fReducedPart=0;}
	203	if(fD0cand) {delete fD0cand; fD0cand=0;}
	204
	205
	206	if(fNofTracks) fNofTracks = 0;
	207
	208	if(fPhiMin) fPhiMin = 0;
	209	if(fPhiMax) fPhiMax = 0;
	210
	211	if(fPtTrigger) fPtTrigger=0;
	212	if(fPhiTrigger) fPhiTrigger=0;
	213	if(fEtaTrigger) fEtaTrigger=0;
	214
	215	if(fDeltaPhi) fDeltaPhi=0;
	216	if(fDeltaEta) fDeltaEta=0;
	217
	218	if(fk0InvMass) fk0InvMass=0;
	219
	220	}
	221	//_____________________________________________________
	222	Bool_t AliHFCorrelator::DefineEventPool(){
	223	// definition of the Pool Manager for Event Mixing
	224
	225
	226	Int_t MaxNofEvents = fhadcuts->GetMaxNEventsInPool();
	227	Int_t MinNofTracks = fhadcuts->GetMinNTracksInPool();
	228	Int_t NofCentBins = fhadcuts->GetNCentPoolBins();
	229	Double_t * CentBins = fhadcuts->GetCentPoolBins();
	230	Int_t NofZVrtxBins = fhadcuts->GetNZvtxPoolBins();
	231	Double_t *ZVrtxBins = fhadcuts->GetZvtxPoolBins();
	232
	233
	234	fPoolMgr = new AliEventPoolManager(MaxNofEvents, MinNofTracks, NofCentBins, CentBins, NofZVrtxBins, ZVrtxBins);
	235	if(!fPoolMgr) return kFALSE;
	236	return kTRUE;
	237	}
	238	//_____________________________________________________
	239	Bool_t AliHFCorrelator::Initialize(){
	240
	241	// std::cout << "AliHFCorrelator::Initialize"<< std::endl;
	242	// AliInfo("AliHFCorrelator::Initialize") ;
	243	if(!fAODEvent){
	244	AliInfo("No AOD event") ;
	245	return kFALSE;
	246	}
	247	//std::cout << "No AOD event" << std::endl;
	248
	249	AliCentrality *centralityObj = 0;
	250	//Int_t multiplicity = -1;
1b8b4a2c	251	//Double_t MultipOrCent = -1;
c683985a	252
c683985a	253	// initialize the pool for event mixing
1b8b4a2c	254	if(!fUseCentrality){ // pp, pA
4640c275	255	//multiplicity = fAODEvent->GetNumberOfTracks();
1b8b4a2c EB	256	//MultipOrCent = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(fAODEvent,-1.,1.);
1b8b4a2c EB	257	fMultCentr = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(fAODEvent,-1.,1.);
c683985a	258	// MultipOrCent = multiplicity; // convert from Int_t to Double_t
	259	// AliInfo(Form("Multiplicity is %f", MultipOrCent));
	260	}
1b8b4a2c EB	261	if(fUseCentrality){ // PbPb
	262	if(!fDMesonCutObject){
	263
0a918d8d	264	centralityObj = ((AliVAODHeader*)fAODEvent->GetHeader())->GetCentralityP();
1b8b4a2c EB	265	fMultCentr = centralityObj->GetCentralityPercentileUnchecked("V0M");
	266	}
	267	else fMultCentr = fDMesonCutObject->GetCentrality(fAODEvent);
c683985a	268	// AliInfo(Form("Centrality is %f", MultipOrCent));
	269	}
	270
	271	AliAODVertex *vtx = fAODEvent->GetPrimaryVertex();
	272	Double_t zvertex = vtx->GetZ(); // zvertex
	273	Double_t * CentBins = fhadcuts->GetCentPoolBins();
	274	Double_t poolmin=CentBins[0];
	275	Double_t poolmax=CentBins[fhadcuts->GetNCentPoolBins()];
	276
	277
1b8b4a2c EB	278	if(TMath::Abs(zvertex)>=10 \|\| fMultCentr>poolmax \|\| fMultCentr < poolmin) {
	279	if(!fUseCentrality)AliInfo(Form("Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,fMultCentr));
	280	if(fUseCentrality) AliInfo(Form("Event with Zvertex = %.2f cm and centrality = %.1f out of pool bounds, SKIPPING",zvertex,fMultCentr));
c683985a	281
	282	return kFALSE;
	283	}
	284
1b8b4a2c	285	fPool = fPoolMgr->GetEventPool(fMultCentr, zvertex);
c683985a	286
c683985a	287	if (!fPool){
1b8b4a2c	288	AliInfo(Form("No pool found for multiplicity = %f, zVtx = %f cm", fMultCentr, zvertex));
c683985a	289	return kFALSE;
	290	}
	291	//fPool->PrintInfo();
	292	return kTRUE;
	293	}
	294
	295	//_____________________________________________________
	296	Bool_t AliHFCorrelator::ProcessEventPool(){
	297	// analysis on Mixed Events
	298	//cout << "AliHFCorrelator::ProcessEventPool"<< endl;
	299	if(!fmixing) return kFALSE;
	300	if(!fPool->IsReady()) return kFALSE;
	301	if(fPool->GetCurrentNEvents()<fhadcuts->GetMinEventsToMix()) return kFALSE;
	302	// fPool->PrintInfo();
	303	fPoolContent = fPool->GetCurrentNEvents();
	304
	305	return kTRUE;
	306
	307	}
	308
	309	//_____________________________________________________
	310	Bool_t AliHFCorrelator::ProcessAssociatedTracks(Int_t EventLoopIndex, const TObjArray* associatedTracks){
	311	// associatedTracks is not deleted, it should be (if needed) deleted in the user task
	312
	313	if(!fmixing){ // analysis on Single Event
	314	if(fAssociatedTracks){
	315	fAssociatedTracks->Delete();
	316	delete fAssociatedTracks;
	317	}
	318	if(fselect==kHadron \|\| fselect ==kKaon){
	319	fAssociatedTracks = AcceptAndReduceTracks(fAODEvent);
	320	fAssociatedTracks->SetOwner(kTRUE);
	321	}
	322	if(fselect==kKZero) {
	323	fAssociatedTracks = AcceptAndReduceKZero(fAODEvent);
	324	fAssociatedTracks->SetOwner(kTRUE);
	325	}
	326	if(fselect==kElectron && associatedTracks) {
	327	fAssociatedTracks=new TObjArray(*associatedTracks);// Maybe better to call the copy constructor
	328	fAssociatedTracks->SetOwner(kFALSE);
	329	}
	330
	331	}
	332
	333	if(fmixing) { // analysis on Mixed Events
	334
	335
	336	fAssociatedTracks = fPool->GetEvent(EventLoopIndex);
	337
	338
	339
	340
	341	} // end if mixing
	342
	343	if(!fAssociatedTracks) return kFALSE;
	344
	345	fNofTracks = fAssociatedTracks->GetEntriesFast();
	346
	347	return kTRUE;
	348
	349	}
	350	//_____________________________________________________
	351	Bool_t AliHFCorrelator::Correlate(Int_t loopindex){
	352
353	if(loopindex >= fNofTracks) return kFALSE;
354	if(!fAssociatedTracks) return kFALSE;
355
356	fReducedPart = (AliReducedParticle*)fAssociatedTracks->At(loopindex);
357
358
359	fDeltaPhi = SetCorrectPhiRange(fPhiTrigger - fReducedPart->Phi());
360
361	fDeltaEta = fEtaTrigger - fReducedPart->Eta();
362
363	return kTRUE;
364
365	}
366
367	//_____________________________________________________
368	Bool_t AliHFCorrelator::PoolUpdate(const TObjArray* associatedTracks){
369
370	if(!fmixing) return kFALSE;
371	if(!fPool) return kFALSE;
372	if(fmixing) { // update the pool for Event Mixing
373	TObjArray* objArr = NULL;
374	if(fselect==kHadron \|\| fselect==kKaon) objArr = (TObjArray*)AcceptAndReduceTracks(fAODEvent);
375	else if(fselect==kKZero) objArr = (TObjArray*)AcceptAndReduceKZero(fAODEvent);
376	else if(fselect==kElectron && associatedTracks){
377	objArr = new TObjArray(*associatedTracks);
378	}
379	else return kFALSE;
380	if(objArr->GetEntriesFast()>0) fPool->UpdatePool(objArr); // updating the pool only if there are entries in the array
381	}
382
383	return kTRUE;
384
385	}
386
387	//_____________________________________________________
388	Double_t AliHFCorrelator::SetCorrectPhiRange(Double_t phi){
389	Double_t pi = TMath::Pi();
390
391	if(phi<fPhiMin) phi = phi + 2*pi;
392	if(phi>fPhiMax) phi = phi - 2*pi;
393
394	return phi;
395	}
396
397	//_____________________________________________________
398	TObjArray* AliHFCorrelator::AcceptAndReduceTracks(AliAODEvent* inputEvent){
399
400	Double_t weight=1.;
4640c275	401	Int_t nTracks = inputEvent->GetNumberOfTracks();
c683985a	402	AliAODVertex * vtx = inputEvent->GetPrimaryVertex();
	403	Double_t pos[3],cov[6];
	404	vtx->GetXYZ(pos);
	405	vtx->GetCovarianceMatrix(cov);
	406	const AliESDVertex vESD(pos,cov,100.,100);
	407
	408	Double_t Bz = inputEvent->GetMagneticField();
	409
	410
	411	TObjArray* tracksClone = new TObjArray;
	412	tracksClone->SetOwner(kTRUE);
	413
	414	//*******************************************************
	415	// use reconstruction
	416	if(fUseReco){
	417	for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
f15c1f69	418	AliAODTrack* track = dynamic_cast<AliAODTrack*>(inputEvent->GetTrack(iTrack));
c683985a	419	if (!track) continue;
	420	if(!fhadcuts->IsHadronSelected(track,&vESD,Bz)) continue; // apply ESD level selections
	421	if(!fhadcuts->Charge(fDCharge,track)) continue; // apply selection on charge, if required
	422
	423	Double_t pT = track->Pt();
	424
	425	//compute impact parameter
	426	Double_t d0z0[2],covd0z0[3];
	427	Double_t d0=-999999.;
	428	if(fUseImpactParameter) track->PropagateToDCA(vtx,Bz,100,d0z0,covd0z0);
	429	else d0z0[0] = 1. ; // random number - be careful with the cuts you applied
	430
	431	if(fUseImpactParameter==1) d0 = TMath::Abs(d0z0[0]); // use impact parameter
	432	if(fUseImpactParameter==2) { // use impact parameter over resolution
	433	if(TMath::Abs(covd0z0[0])>0.00000001) d0 = TMath::Abs(d0z0[0])/TMath::Sqrt(covd0z0[0]);
	434	else d0 = -1.; // if the resoultion is Zero, rejects the track - to be on the safe side
	435
	436	}
	437
	438	if(fmontecarlo) {// THIS TO BE CHECKED
	439	Int_t hadLabel = track->GetLabel();
	440	if(hadLabel < 0) continue;
	441	}
	442
	443	if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
	444	Bool_t rejectsoftpi = kTRUE;// TO BE CHECKED: DO WE WANT IT TO kTRUE AS A DEFAULT?
	445	if(fD0cand && !fmixing) rejectsoftpi = fhadcuts->InvMassDstarRejection(fD0cand,track,fhypD0); // TO BE CHECKED: WHY NOT FOR EM?
	446
	447
	448	if(fselect ==kKaon){
	449	if(!fhadcuts->CheckKaonCompatibility(track,fmontecarlo,fmcArray,fPIDmode)) continue; // check if it is a Kaon - data and MC
	450	}
	451	weight=fhadcuts->GetTrackWeight(pT,track->Eta(),pos[2]);
	452	tracksClone->Add(new AliReducedParticle(track->Eta(), track->Phi(), pT,track->GetLabel(),track->GetID(),d0,rejectsoftpi,track->Charge(),weight));
	453	} // end loop on tracks
	454	} // end if use reconstruction kTRUE
	455
	456	//*******************************************************
	457
	458	//use MC truth
	459	if(fmontecarlo && !fUseReco){
	460
	461	for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) {
	462	AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));
	463	if (!mcPart) {
	464	AliWarning("MC Particle not found in tree, skipping");
	465	continue;
	466	}
	467	if(!mcPart->Charge()) continue; // consider only charged tracks
	468
	469	Int_t PDG =TMath::Abs(mcPart->PdgCode());
	470	if(fselect ==kHadron) {if(!((PDG==321)\|\|(PDG==211)\|\|(PDG==2212)\|\|(PDG==13)\|\|(PDG==11))) continue;} // select only if kaon, pion, proton, muon or electron
	471	else if(fselect ==kKaon) {if(!(PDG==321)) continue;} // select only if kaon
	472	else if(fselect ==kElectron) {if(!(PDG==11)) continue;} // select only if electron
	473
	474	Double_t pT = mcPart->Pt();
	475	Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet
	476	if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
	477
	478	tracksClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,iPart,-1,d0,kFALSE,mcPart->Charge()));
	479	}
	480
	481	} // end if use MC truth
	482
483
484	return tracksClone;
485	}
486
487	//_____________________________________________________
488	TObjArray* AliHFCorrelator::AcceptAndReduceKZero(AliAODEvent* inputEvent){
489
490	Int_t nOfVZeros = inputEvent->GetNumberOfV0s();
491	TObjArray* KZeroClone = new TObjArray;
492	AliAODVertex vertex1 = (AliAODVertex)inputEvent->GetPrimaryVertex();
493
494	// use reconstruction
495	if(fUseReco){
496	Int_t v0label = -1;
497	Int_t pdgDgK0toPipi[2] = {211,211};
498	Double_t mPDGK0=0.497614;//TDatabasePDG::Instance()->GetParticle(310)->Mass();
499	const Int_t size = inputEvent->GetNumberOfV0s();
500	Int_t prevnegID[size];
501	Int_t prevposID[size];
502	for(Int_t iv0 =0; iv0< nOfVZeros; iv0++){// loop on all v0 candidates
503	AliAODv0 v0 = (static_cast<AliAODEvent> (inputEvent))->GetV0(iv0);
504	if(!v0) {
505	AliInfo(Form("is not a v0 at step, %d", iv0)) ;
506	//cout << "is not a v0 at step " << iv0 << endl;
507	continue;
508	}
509
510	if(!fhadcuts->IsKZeroSelected(v0,vertex1)) continue; // check if it is a k0
511
512	// checks to avoid double counting
513	Int_t negID = -999;
514	Int_t posID = -998;
515	//Int_t a = 0;
516	prevnegID[iv0] = -997;
517	prevposID[iv0]= -996;
518	negID = v0->GetNegID();
519	posID = v0->GetPosID();
520	Bool_t DoubleCounting = kFALSE;
521
522	for(Int_t k=0; k<iv0; k++){
523	if(((negID==prevnegID[k])&&(posID==prevposID[k]))\|\|((negID==prevposID[k])&&(posID==prevnegID[k]))){
524	DoubleCounting = kTRUE;
525	//a=k;
526	break;
527	}//end if
528	} // end for
529
530	if(DoubleCounting) continue;
531	else{
532	prevposID[iv0] = posID;
533	prevnegID[iv0] = negID;
534	}
535
536	if(fmontecarlo) v0label = v0->MatchToMC(310,fmcArray, 0, pdgDgK0toPipi); //match a K0 short
537	Double_t k0pt = v0->Pt();
538	Double_t k0eta = v0->Eta();
539	Double_t k0Phi = v0->Phi();
540	fk0InvMass = v0->MassK0Short();
541
542	//if (loopindex == 0) {
543	// if(!plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectra"))->Fill(k0InvMass,k0pt); // spectra for all k0
544	// if(plotassociation) ((TH2F)fOutput->FindObject("KZeroSpectraifHF"))->Fill(k0InvMass,k0pt); // spectra for k0 in association with a D
545	//}
546	// if there are more D* candidates per event, loopindex == 0 makes sure you fill the mass spectra only once!
547
548	if(TMath::Abs(fk0InvMass-mPDGK0)>3*0.004) continue; // select candidates within 3 sigma
549	KZeroClone->Add(new AliReducedParticle(k0eta,k0Phi,k0pt,v0label));
550
551	}
552	} // end if use reconstruction kTRUE
553
554
555
556	//*********************************************************************//
557	//use MC truth
558	if(fmontecarlo && !fUseReco){
559
560	for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) {
561	AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));
562	if (!mcPart) {
563	AliWarning("MC Particle not found in tree, skipping");
564	continue;
565	}
566
567	Int_t PDG =TMath::Abs(mcPart->PdgCode());
568	if(!(PDG==310)) continue; // select only if k0 short
569
570	Double_t pT = mcPart->Pt();
571	Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet
572	if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
573
574	KZeroClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,iPart,-1,d0,kTRUE,mcPart->Charge()));
575	}
576
577	} // end if use MC truth
578
579
580
581	return KZeroClone;
582	}