[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisUtils.cxx

#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliVVertex.h"
#include "AliLog.h"
#include "AliAODVertex.h"
#include "AliVTrack.h"
#include "AliVEvent.h"
#include <TMatrixDSym.h>
#include <TMath.h>

#include "AliAnalysisUtils.h"

ClassImp(AliAnalysisUtils)

//______________________________________________________________________
AliAnalysisUtils::AliAnalysisUtils():TObject(),
  fisAOD(kTRUE),
  fMinVtxContr(0),
  fMaxVtxZ(10.),
  fCutOnZVertexSPD(kTRUE),
  fUseMVPlpSelection(kFALSE),
  fUseOutOfBunchPileUp(kFALSE),
  fMinPlpContribMV(5),
  fMaxPlpChi2MV(5.),
  fMinWDistMV(15.),
  fCheckPlpFromDifferentBCMV(kFALSE),
  fMinPlpContribSPD(5),
  fMinPlpZdistSPD(0.8),
  fnSigmaPlpZdistSPD(3.),
  fnSigmaPlpDiamXYSPD(2.),
  fnSigmaPlpDiamZSPD(5.),
  fUseSPDCutInMultBins(kFALSE)
{
  // Default contructor
}

//______________________________________________________________________
Bool_t AliAnalysisUtils::IsVertexSelected2013pA(AliVEvent *event)
{
  Bool_t accept = kFALSE;
  
  // Check whether the event is of AOD or ESD type
  AliAODEvent *aod = 0x0;
  AliESDEvent *esd = 0x0;
  aod = dynamic_cast<AliAODEvent*>(event);
  esd = dynamic_cast<AliESDEvent*>(event);

  if(aod) { 
    fisAOD = kTRUE; 
  } else {
    fisAOD = kFALSE;
    if(!esd) {
      AliFatal("Event is neither of AOD nor ESD type");
      return accept;
    }
  }

  const AliVVertex *trkVtx = fisAOD ? 
    dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) : 
    dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertex()) ;
  if(!trkVtx || trkVtx->GetNContributors()<=fMinVtxContr){
    accept = kFALSE;
    return accept;
  }

  TString vtxTtl = trkVtx->GetTitle();
  if (!vtxTtl.Contains("VertexerTracks")) return accept;

  Float_t zvtx = trkVtx->GetZ();
  const AliVVertex* spdVtx = fisAOD ? 
    dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertexSPD()) : 
    dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexSPD()) ;
  if (spdVtx->GetNContributors()<=fMinVtxContr) return accept;

  TString vtxTyp = spdVtx->GetTitle();
  Double_t cov[6]={0};
  spdVtx->GetCovarianceMatrix(cov);
  Double_t zRes = TMath::Sqrt(cov[5]);
  if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return accept;
  if (fCutOnZVertexSPD && TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return accept;

  if (TMath::Abs(zvtx) > fMaxVtxZ) return accept;

  return kTRUE;
}

//______________________________________________________________________
Bool_t AliAnalysisUtils::IsFirstEventInChunk(AliVEvent *event)
{

  Bool_t accept = kFALSE;

  // Check whether the event is of AOD or ESD type
  AliAODEvent *aod = 0x0;
  AliESDEvent *esd = 0x0;
  aod = dynamic_cast<AliAODEvent*>(event);
  esd = dynamic_cast<AliESDEvent*>(event);

  if(aod) { 
    fisAOD = kTRUE; 
  } else {
    fisAOD = kFALSE;
    if(!esd) {
      AliFatal("Event is neither of AOD nor ESD type");
      return accept;
    }
  }

  if(fisAOD){
    AliAODHeader *aodheader = 0x0;
    aodheader = aod->GetHeader();
    if(!aodheader){
      AliFatal("AOD header not there ?!");
      return kFALSE;
    }
    if(aodheader->GetEventNumberESDFile()==0) accept = kTRUE;
  } else {
    if(esd->GetEventNumberInFile()==0) accept = kTRUE;
  }
  
  return accept;
}

//______________________________________________________________________
Bool_t AliAnalysisUtils::IsPileUpEvent(AliVEvent *event)
{
  Bool_t isPileUp=kFALSE;
  //check for multiple vertices
  if(fUseMVPlpSelection)isPileUp=IsPileUpMV(event);
  else isPileUp=IsPileUpSPD(event);
  //check for different BC 
  if(fUseOutOfBunchPileUp && IsOutOfBunchPileUp(event))isPileUp=kTRUE;
  
  return isPileUp;
}

//______________________________________________________________________
Bool_t AliAnalysisUtils::IsPileUpMV(AliVEvent *event)
{
  // check for multi-vertexer pile-up
  const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
  const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
  //
  if (!aod && !esd) {
    AliFatal("Event is neither of AOD nor ESD type");
    return kFALSE;
  }
  //
  const AliVVertex* vtPrm = 0;
  const AliVVertex* vtPlp = 0;
  Int_t nPlp = 0;
  //
  if (aod) {
    if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
    vtPrm = aod->GetPrimaryVertex();
    if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
  }
  else {
    if ( !(nPlp=esd->GetNumberOfPileupVerticesTracks())) return kFALSE;
    vtPrm = esd->GetPrimaryVertexTracks();
    if (((AliESDVertex*)vtPrm)->GetStatus()!=1) return kTRUE; // there are pile-up vertices but no primary
  }
  Int_t bcPrim = vtPrm->GetBC();
  //
  for (Int_t ipl=0;ipl<nPlp;ipl++) {
    vtPlp = aod ? (const AliVVertex*)aod->GetPileupVertexTracks(ipl) : (const AliVVertex*)esd->GetPileupVertexTracks(ipl);
    //
    if (vtPlp->GetNContributors() < fMinPlpContribMV) continue;
    if (vtPlp->GetChi2perNDF() > fMaxPlpChi2MV) continue;
    if(fCheckPlpFromDifferentBCMV)
      {
	Int_t bcPlp = vtPlp->GetBC();
	if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
      }
    //
    Double_t wDst = GetWDist(vtPrm,vtPlp);
    if (wDst<fMinWDistMV) continue;
    //
    return kTRUE; // pile-up: well separated vertices
  }
  //
  return kFALSE;
  //
}

//______________________________________________________________________
Bool_t AliAnalysisUtils::IsPileUpSPD(AliVEvent *event)
{
  // check for SPD pile-up
  const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
  const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
  //
  if (!aod && !esd) {
    AliFatal("Event is neither of AOD nor ESD type");
    return kFALSE;
  }
  //
  if (aod) return (fUseSPDCutInMultBins)?aod->IsPileupFromSPDInMultBins():aod->IsPileupFromSPD(fMinPlpContribSPD,fMinPlpZdistSPD,fnSigmaPlpZdistSPD,fnSigmaPlpDiamXYSPD,fnSigmaPlpDiamZSPD);
  else return (fUseSPDCutInMultBins)?esd->IsPileupFromSPDInMultBins():esd->IsPileupFromSPD(fMinPlpContribSPD,fMinPlpZdistSPD,fnSigmaPlpZdistSPD,fnSigmaPlpDiamXYSPD,fnSigmaPlpDiamZSPD);
}

//______________________________________________________________________
Bool_t AliAnalysisUtils::IsOutOfBunchPileUp(AliVEvent *event)
{
  // check for SPD pile-up
  const AliAODEvent *aod = dynamic_cast<const AliAODEvent*>(event);
  const AliESDEvent *esd = dynamic_cast<const AliESDEvent*>(event);
  //
  if (!aod && !esd) {
    AliFatal("Event is neither of AOD nor ESD type");
    return kFALSE;
  }
  Int_t bc2 = (aod)?aod->GetHeader()->GetIRInt2ClosestInteractionMap():esd->GetHeader()->GetIRInt2ClosestInteractionMap();
  if (bc2 != 0)
    return kTRUE;
  
  Int_t bc1 = (aod)?aod->GetHeader()->GetIRInt1ClosestInteractionMap():esd->GetHeader()->GetIRInt1ClosestInteractionMap();
  if (bc1 != 0)
    return kTRUE;
  
  return kFALSE;
}

//______________________________________________________________________
Double_t AliAnalysisUtils::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
{
  // calculate sqrt of weighted distance to other vertex
  if (!v0 || !v1) {
    printf("One of vertices is not valid\n");
    return 0;
  }
  static TMatrixDSym vVb(3);
  Double_t dist = -1;
  Double_t dx = v0->GetX()-v1->GetX();
  Double_t dy = v0->GetY()-v1->GetY();
  Double_t dz = v0->GetZ()-v1->GetZ();
  Double_t cov0[6],cov1[6];
  v0->GetCovarianceMatrix(cov0);
  v1->GetCovarianceMatrix(cov1);
  vVb(0,0) = cov0[0]+cov1[0];
  vVb(1,1) = cov0[2]+cov1[2];
  vVb(2,2) = cov0[5]+cov1[5];
  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
  vVb.InvertFast();
  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
    +    2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
  return dist>0 ? TMath::Sqrt(dist) : -1; 

}
Commit	Line	Data
77c42b43	1	#include "AliVEvent.h"
	2	#include "AliESDEvent.h"
	3	#include "AliAODEvent.h"
	4	#include "AliVVertex.h"
	5	#include "AliLog.h"
	6	#include "AliAODVertex.h"
c9df49b3	7	#include "AliVTrack.h"
	8	#include "AliVEvent.h"
	9	#include <TMatrixDSym.h>
	10	#include <TMath.h>
77c42b43	11
	12	#include "AliAnalysisUtils.h"
	13
	14	ClassImp(AliAnalysisUtils)
	15
	16	//______________________________________________________________________
	17	AliAnalysisUtils::AliAnalysisUtils():TObject(),
	18	fisAOD(kTRUE),
	19	fMinVtxContr(0),
	20	fMaxVtxZ(10.),
c9df49b3	21	fCutOnZVertexSPD(kTRUE),
	22	fUseMVPlpSelection(kFALSE),
	23	fUseOutOfBunchPileUp(kFALSE),
	24	fMinPlpContribMV(5),
	25	fMaxPlpChi2MV(5.),
	26	fMinWDistMV(15.),
	27	fCheckPlpFromDifferentBCMV(kFALSE),
	28	fMinPlpContribSPD(5),
	29	fMinPlpZdistSPD(0.8),
	30	fnSigmaPlpZdistSPD(3.),
	31	fnSigmaPlpDiamXYSPD(2.),
	32	fnSigmaPlpDiamZSPD(5.),
	33	fUseSPDCutInMultBins(kFALSE)
77c42b43	34	{
	35	// Default contructor
	36	}
	37
	38	//______________________________________________________________________
	39	Bool_t AliAnalysisUtils::IsVertexSelected2013pA(AliVEvent *event)
	40	{
	41	Bool_t accept = kFALSE;
c9df49b3	42
77c42b43	43	// Check whether the event is of AOD or ESD type
	44	AliAODEvent *aod = 0x0;
	45	AliESDEvent *esd = 0x0;
	46	aod = dynamic_cast<AliAODEvent*>(event);
	47	esd = dynamic_cast<AliESDEvent*>(event);
	48
	49	if(aod) {
	50	fisAOD = kTRUE;
	51	} else {
	52	fisAOD = kFALSE;
	53	if(!esd) {
	54	AliFatal("Event is neither of AOD nor ESD type");
	55	return accept;
	56	}
	57	}
	58
	59	const AliVVertex *trkVtx = fisAOD ?
	60	dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) :
	61	dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertex()) ;
	62	if(!trkVtx \|\| trkVtx->GetNContributors()<=fMinVtxContr){
	63	accept = kFALSE;
	64	return accept;
	65	}
	66
	67	TString vtxTtl = trkVtx->GetTitle();
	68	if (!vtxTtl.Contains("VertexerTracks")) return accept;
	69
	70	Float_t zvtx = trkVtx->GetZ();
	71	const AliVVertex* spdVtx = fisAOD ?
	72	dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertexSPD()) :
	73	dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexSPD()) ;
	74	if (spdVtx->GetNContributors()<=fMinVtxContr) return accept;
	75
	76	TString vtxTyp = spdVtx->GetTitle();
	77	Double_t cov[6]={0};
	78	spdVtx->GetCovarianceMatrix(cov);
	79	Double_t zRes = TMath::Sqrt(cov[5]);
	80	if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return accept;
	81	if (fCutOnZVertexSPD && TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return accept;
	82
	83	if (TMath::Abs(zvtx) > fMaxVtxZ) return accept;
	84
	85	return kTRUE;
	86	}
	87
	88	//______________________________________________________________________
	89	Bool_t AliAnalysisUtils::IsFirstEventInChunk(AliVEvent *event)
	90	{
	91
	92	Bool_t accept = kFALSE;
	93
	94	// Check whether the event is of AOD or ESD type
	95	AliAODEvent *aod = 0x0;
	96	AliESDEvent *esd = 0x0;
	97	aod = dynamic_cast<AliAODEvent*>(event);
	98	esd = dynamic_cast<AliESDEvent*>(event);
	99
	100	if(aod) {
	101	fisAOD = kTRUE;
	102	} else {
	103	fisAOD = kFALSE;
	104	if(!esd) {
	105	AliFatal("Event is neither of AOD nor ESD type");
	106	return accept;
107	}
108	}
109
110	if(fisAOD){
111	AliAODHeader *aodheader = 0x0;
112	aodheader = aod->GetHeader();
113	if(!aodheader){
114	AliFatal("AOD header not there ?!");
115	return kFALSE;
116	}
117	if(aodheader->GetEventNumberESDFile()==0) accept = kTRUE;
118	} else {
119	if(esd->GetEventNumberInFile()==0) accept = kTRUE;
120	}
c9df49b3	121
77c42b43	122	return accept;
77c42b43	123	}
c9df49b3	124
	125	//______________________________________________________________________
	126	Bool_t AliAnalysisUtils::IsPileUpEvent(AliVEvent *event)
	127	{
	128	Bool_t isPileUp=kFALSE;
	129	//check for multiple vertices
	130	if(fUseMVPlpSelection)isPileUp=IsPileUpMV(event);
	131	else isPileUp=IsPileUpSPD(event);
	132	//check for different BC
	133	if(fUseOutOfBunchPileUp && IsOutOfBunchPileUp(event))isPileUp=kTRUE;
	134
	135	return isPileUp;
	136	}
	137
	138	//______________________________________________________________________
	139	Bool_t AliAnalysisUtils::IsPileUpMV(AliVEvent *event)
	140	{
	141	// check for multi-vertexer pile-up
	142	const AliAODEvent aod = dynamic_cast<const AliAODEvent>(event);
	143	const AliESDEvent esd = dynamic_cast<const AliESDEvent>(event);
	144	//
	145	if (!aod && !esd) {
	146	AliFatal("Event is neither of AOD nor ESD type");
	147	return kFALSE;
	148	}
	149	//
	150	const AliVVertex* vtPrm = 0;
	151	const AliVVertex* vtPlp = 0;
	152	Int_t nPlp = 0;
	153	//
	154	if (aod) {
	155	if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
	156	vtPrm = aod->GetPrimaryVertex();
	157	if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
	158	}
	159	else {
	160	if ( !(nPlp=esd->GetNumberOfPileupVerticesTracks())) return kFALSE;
	161	vtPrm = esd->GetPrimaryVertexTracks();
	162	if (((AliESDVertex*)vtPrm)->GetStatus()!=1) return kTRUE; // there are pile-up vertices but no primary
	163	}
	164	Int_t bcPrim = vtPrm->GetBC();
	165	//
	166	for (Int_t ipl=0;ipl<nPlp;ipl++) {
	167	vtPlp = aod ? (const AliVVertex)aod->GetPileupVertexTracks(ipl) : (const AliVVertex)esd->GetPileupVertexTracks(ipl);
	168	//
	169	if (vtPlp->GetNContributors() < fMinPlpContribMV) continue;
	170	if (vtPlp->GetChi2perNDF() > fMaxPlpChi2MV) continue;
	171	if(fCheckPlpFromDifferentBCMV)
	172	{
	173	Int_t bcPlp = vtPlp->GetBC();
	174	if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
	175	}
	176	//
	177	Double_t wDst = GetWDist(vtPrm,vtPlp);
	178	if (wDst<fMinWDistMV) continue;
	179	//
	180	return kTRUE; // pile-up: well separated vertices
	181	}
	182	//
	183	return kFALSE;
	184	//
	185	}
	186
	187	//______________________________________________________________________
188	Bool_t AliAnalysisUtils::IsPileUpSPD(AliVEvent *event)
189	{
190	// check for SPD pile-up
191	const AliAODEvent aod = dynamic_cast<const AliAODEvent>(event);
192	const AliESDEvent esd = dynamic_cast<const AliESDEvent>(event);
193	//
194	if (!aod && !esd) {
195	AliFatal("Event is neither of AOD nor ESD type");
196	return kFALSE;
197	}
198	//
199	if (aod) return (fUseSPDCutInMultBins)?aod->IsPileupFromSPDInMultBins():aod->IsPileupFromSPD(fMinPlpContribSPD,fMinPlpZdistSPD,fnSigmaPlpZdistSPD,fnSigmaPlpDiamXYSPD,fnSigmaPlpDiamZSPD);
200	else return (fUseSPDCutInMultBins)?esd->IsPileupFromSPDInMultBins():esd->IsPileupFromSPD(fMinPlpContribSPD,fMinPlpZdistSPD,fnSigmaPlpZdistSPD,fnSigmaPlpDiamXYSPD,fnSigmaPlpDiamZSPD);
201	}
202
203	//______________________________________________________________________
204	Bool_t AliAnalysisUtils::IsOutOfBunchPileUp(AliVEvent *event)
205	{
206	// check for SPD pile-up
207	const AliAODEvent aod = dynamic_cast<const AliAODEvent>(event);
208	const AliESDEvent esd = dynamic_cast<const AliESDEvent>(event);
209	//
210	if (!aod && !esd) {
211	AliFatal("Event is neither of AOD nor ESD type");
212	return kFALSE;
213	}
214	Int_t bc2 = (aod)?aod->GetHeader()->GetIRInt2ClosestInteractionMap():esd->GetHeader()->GetIRInt2ClosestInteractionMap();
215	if (bc2 != 0)
216	return kTRUE;
217
218	Int_t bc1 = (aod)?aod->GetHeader()->GetIRInt1ClosestInteractionMap():esd->GetHeader()->GetIRInt1ClosestInteractionMap();
219	if (bc1 != 0)
220	return kTRUE;
221
222	return kFALSE;
223	}
224
225	//______________________________________________________________________
226	Double_t AliAnalysisUtils::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
227	{
228	// calculate sqrt of weighted distance to other vertex
229	if (!v0 \|\| !v1) {
230	printf("One of vertices is not valid\n");
231	return 0;
232	}
233	static TMatrixDSym vVb(3);
234	Double_t dist = -1;
235	Double_t dx = v0->GetX()-v1->GetX();
236	Double_t dy = v0->GetY()-v1->GetY();
237	Double_t dz = v0->GetZ()-v1->GetZ();
238	Double_t cov0[6],cov1[6];
239	v0->GetCovarianceMatrix(cov0);
240	v1->GetCovarianceMatrix(cov1);
241	vVb(0,0) = cov0[0]+cov1[0];
242	vVb(1,1) = cov0[2]+cov1[2];
243	vVb(2,2) = cov0[5]+cov1[5];
244	vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
245	vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
246	vVb.InvertFast();
247	if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
248	dist = vVb(0,0)dxdx + vVb(1,1)dydy + vVb(2,2)dzdz
249	+ 2vVb(0,1)dxdy + 2vVb(0,2)dxdz + 2vVb(1,2)dy*dz;
250	return dist>0 ? TMath::Sqrt(dist) : -1;
251
252	}