[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAODRecoDecayHF.cxx

/**************************************************************************
 * Copyright(c) 1998-2006, 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.                  *
 **************************************************************************/

/* $Id$ */

/////////////////////////////////////////////////////////////
//
// Base class for AOD reconstructed heavy-flavour decay
//
// Author: A.Dainese, andrea.dainese@lnl.infn.it
/////////////////////////////////////////////////////////////

#include <TDatabasePDG.h>
#include <TVector3.h>
#include <TRandom.h>
#include "AliAODRecoDecay.h"
#include "AliAODRecoDecayHF.h"
#include "AliAODEvent.h"
#include "AliVertexerTracks.h"
#include "AliExternalTrackParam.h"
#include "AliKFVertex.h"
#include "AliVVertex.h"
#include "AliESDVertex.h"

ClassImp(AliAODRecoDecayHF)

//--------------------------------------------------------------------------
AliAODRecoDecayHF::AliAODRecoDecayHF() :
  AliAODRecoDecay(),
  fOwnPrimaryVtx(0x0),
  fEventPrimaryVtx(),
  fListOfCuts(),
  fd0err(0x0), 
  fProngID(0x0),
  fSelectionMap(0)
{
  //
  // Default Constructor
  //
}
//--------------------------------------------------------------------------
AliAODRecoDecayHF::AliAODRecoDecayHF(AliAODVertex *vtx2,Int_t nprongs,Short_t charge,
				     Double_t *px,Double_t *py,Double_t *pz,
				     Double_t *d0,Double_t *d0err) :
  AliAODRecoDecay(vtx2,nprongs,charge,px,py,pz,d0),
  fOwnPrimaryVtx(0x0),
  fEventPrimaryVtx(),
  fListOfCuts(),
  fd0err(0x0),
  fProngID(0x0),
  fSelectionMap(0)
{
  //
  // Constructor with AliAODVertex for decay vertex
  //
  fd0err = new Double_t[GetNProngs()];
  for(Int_t i=0; i<GetNProngs(); i++) fd0err[i] = d0err[i];
}
//--------------------------------------------------------------------------
AliAODRecoDecayHF::AliAODRecoDecayHF(AliAODVertex *vtx2,Int_t nprongs,Short_t charge,
				     Double_t *d0,Double_t *d0err) :
  AliAODRecoDecay(vtx2,nprongs,charge,d0),
  fOwnPrimaryVtx(0x0),
  fEventPrimaryVtx(),
  fListOfCuts(),
  fd0err(0x0),
  fProngID(0x0),
  fSelectionMap(0)
{
  //
  // Constructor with AliAODVertex for decay vertex and without prongs momenta
  //
  fd0err = new Double_t[GetNProngs()];
  for(Int_t i=0; i<GetNProngs(); i++) fd0err[i] = d0err[i];
}
//--------------------------------------------------------------------------
AliAODRecoDecayHF::AliAODRecoDecayHF(Double_t vtx1[3],Double_t vtx2[3],
				     Int_t nprongs,Short_t charge,
				     Double_t *px,Double_t *py,Double_t *pz,
				     Double_t *d0) :
  AliAODRecoDecay(0x0,nprongs,charge,px,py,pz,d0),
  fOwnPrimaryVtx(0x0),
  fEventPrimaryVtx(),
  fListOfCuts(),
  fd0err(0x0),
  fProngID(0x0), 
  fSelectionMap(0)
{
  //
  // Constructor that can used for a "MC" object
  //

  fOwnPrimaryVtx = new AliAODVertex(vtx1);

  AliAODVertex *vtx = new AliAODVertex(vtx2);
  SetOwnSecondaryVtx(vtx);

}
//--------------------------------------------------------------------------
AliAODRecoDecayHF::AliAODRecoDecayHF(const AliAODRecoDecayHF &source) :
  AliAODRecoDecay(source),
  fOwnPrimaryVtx(0x0),
  fEventPrimaryVtx(source.fEventPrimaryVtx),
  fListOfCuts(source.fListOfCuts),
  fd0err(0x0),
  fProngID(0x0),
  fSelectionMap(source.fSelectionMap)
{
  //
  // Copy constructor
  //
  if(source.GetOwnPrimaryVtx()) fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));

  if(source.GetNProngs()>0) {
    fd0err = new Double_t[GetNProngs()];
    memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
    if(source.fProngID) {
      fProngID = new UShort_t[GetNProngs()];
      memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
    }
  }
}
//--------------------------------------------------------------------------
AliAODRecoDecayHF &AliAODRecoDecayHF::operator=(const AliAODRecoDecayHF &source)
{
  //
  // assignment operator
  //
  if(&source == this) return *this;

  AliAODRecoDecay::operator=(source);

  fEventPrimaryVtx = source.fEventPrimaryVtx;
  fListOfCuts = source.fListOfCuts;
  fSelectionMap = source.fSelectionMap;

  if(source.GetOwnPrimaryVtx()) {
    delete fOwnPrimaryVtx;
    fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
  }

  if(source.GetNProngs()>0) {
    if(source.fd0err) {
      delete [] fd0err;
      fd0err = new Double_t[GetNProngs()];
      memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
    }
    if(source.fProngID) {
      delete [] fProngID;
      fProngID = new UShort_t[GetNProngs()];
      memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
    }
  }
  return *this;
}
//--------------------------------------------------------------------------
AliAODRecoDecayHF::~AliAODRecoDecayHF() {
  //  
  // Default Destructor
  //
  if(fOwnPrimaryVtx) delete fOwnPrimaryVtx;
  if(fd0err) delete [] fd0err;
  if(fProngID) delete [] fProngID;
}
//---------------------------------------------------------------------------
AliKFParticle *AliAODRecoDecayHF::ApplyVertexingKF(Int_t *iprongs,Int_t nprongs,Int_t *pdgs,Bool_t topoCostraint, Double_t bzkG, Double_t *mass) const {
  //
  // Applies the KF vertexer 
  // Int_t iprongs[nprongs] = indices of the prongs to be used from the vertexer
  // Int_t pdgs[nprongs] = pdgs assigned to the prongs, needed to define the AliKFParticle
  // Bool_t topoCostraint = if kTRUE, the topological constraint is applied
  // Double_t bzkG = magnetic field
  // Double_t mass[2] = {mass, sigma} for the mass constraint (if mass[0]>0 the constraint is applied).
  //

  AliKFParticle::SetField(bzkG);
  AliKFParticle *vertexKF=0;
  
  AliKFVertex copyKF;
  Int_t nt=0,ntcheck=0;

  Double_t pos[3]={0.,0.,0.};
  if(!fOwnPrimaryVtx) {
    printf("AliAODRecoDecayHF::ApplyVertexingKF(): cannot apply because primary vertex is not found\n");
    return vertexKF;
  }
  fOwnPrimaryVtx->GetXYZ(pos);
  Int_t contr=fOwnPrimaryVtx->GetNContributors();
  Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
  fOwnPrimaryVtx->GetCovarianceMatrix(covmatrix);
  Double_t chi2=fOwnPrimaryVtx->GetChi2();
  AliESDVertex primaryVtx2(pos,covmatrix,chi2,contr,"Vertex");
 

  if(topoCostraint){
   copyKF=AliKFVertex(primaryVtx2);
   nt=primaryVtx2.GetNContributors();
   ntcheck=nt;
  }

  vertexKF = new AliKFParticle();
  for(Int_t i= 0;i<nprongs;i++){
    Int_t ipr=iprongs[i];
    AliAODTrack *aodTrack = (AliAODTrack*)GetDaughter(ipr);
    if(!aodTrack) {
      printf("AliAODRecoDecayHF::ApplyVertexingKF(): no daughters available\n");
      delete vertexKF; vertexKF=NULL;
      return vertexKF;
    }
    AliKFParticle daughterKF(*aodTrack,pdgs[i]);
    vertexKF->AddDaughter(daughterKF);
    
    if(topoCostraint && nt>0){
      //Int_t index=(Int_t)GetProngID(ipr);
      if(!aodTrack->GetUsedForPrimVtxFit()) continue;
      copyKF -= daughterKF;
      ntcheck--;
    }
  }
  
  if(topoCostraint){
    if(ntcheck>0) {
      copyKF += (*vertexKF);
      vertexKF->SetProductionVertex(copyKF);
    }
 }
  
  if(mass[0]>0.){
    vertexKF->SetMassConstraint(mass[0],mass[1]);
  }
  
  return vertexKF;
}
//---------------------------------------------------------------------------
AliAODVertex* AliAODRecoDecayHF::RemoveDaughtersFromPrimaryVtx(AliAODEvent *aod) {
  //
  // This method returns a primary vertex without the daughter tracks of the 
  // candidate and it recalculates the impact parameters and errors.
  // 
  // The output vertex is created with "new". The user has to 
  // set it to the candidate with SetOwnPrimaryVtx(), unset it at the end 
  // of processing with UnsetOwnPrimaryVtx() and delete it.
  // If a NULL pointer is returned, the removal failed (too few tracks left).
  //
  // For the moment, the primary vertex is recalculated from scratch without
  // the daughter tracks.
  //

  AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
  if(!vtxAOD) return 0;
  TString title=vtxAOD->GetTitle();
  if(!title.Contains("VertexerTracks")) return 0;


  AliVertexerTracks *vertexer = new AliVertexerTracks(aod->GetMagneticField());

  Int_t ndg = GetNDaughters();

  vertexer->SetITSMode();
  vertexer->SetMinClusters(3);
  vertexer->SetConstraintOff();

  if(title.Contains("WithConstraint")) {
    Float_t diamondcovxy[3];
    aod->GetDiamondCovXY(diamondcovxy);
    Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
    vertexer->SetVtxStart(diamond);
    delete diamond; diamond=NULL;
  }

  Int_t skipped[10]; for(Int_t i=0;i<10;i++) skipped[i]=-1;
  Int_t nTrksToSkip=0,id;
  AliAODTrack *t = 0;
  for(Int_t i=0; i<ndg; i++) {
    t = (AliAODTrack*)GetDaughter(i);
    id = (Int_t)t->GetID();
    if(id<0) continue;
    skipped[nTrksToSkip++] = id;
  }
  vertexer->SetSkipTracks(nTrksToSkip,skipped);
  AliESDVertex *vtxESDNew = vertexer->FindPrimaryVertex(aod);

  delete vertexer; vertexer=NULL;

  if(!vtxESDNew) return 0;
  if(vtxESDNew->GetNContributors()<=0) { 
    delete vtxESDNew; vtxESDNew=NULL;
    return 0;
  }

  // convert to AliAODVertex
  Double_t pos[3],cov[6],chi2perNDF;
  vtxESDNew->GetXYZ(pos); // position
  vtxESDNew->GetCovMatrix(cov); //covariance matrix
  chi2perNDF = vtxESDNew->GetChi2toNDF();
  delete vtxESDNew; vtxESDNew=NULL;

  AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);

  RecalculateImpPars(vtxAODNew,aod);

  return vtxAODNew;
}
//-----------------------------------------------------------------------------------
void AliAODRecoDecayHF::RecalculateImpPars(AliAODVertex *vtxAODNew,AliAODEvent* aod) {
  //
  // now recalculate the daughters impact parameters
  //
  Double_t dz[2],covdz[3];
  for(Int_t i=0; i<GetNDaughters(); i++) {
    AliAODTrack *t = (AliAODTrack*)GetDaughter(i);
    AliExternalTrackParam etp; etp.CopyFromVTrack(t);
    if(etp.PropagateToDCA(vtxAODNew,aod->GetMagneticField(),3.,dz,covdz)) {
      fd0[i]    = dz[0];
      fd0err[i] = TMath::Sqrt(covdz[0]);
    }
  }

  return;
}
//-----------------------------------------------------------------------------------
void AliAODRecoDecayHF::Misalign(TString misal) {
  //
  // Method to smear the impact parameter of the duaghter tracks
  // and the sec. vtx position accordinlgy 
  // Useful to study effect of misalignment.
  // The starting point are parameterisations of the impact parameter resolution
  // from MC and data 
  // Errors on d0 and vtx are not recalculated (to be done)
  //
  if(misal=="null")return;
  Double_t pard0rphiMC[3]={36.7,36.,1.25};// d0(pt)=[0]+[1]/(pt^[2]); in micron, conversion to cm is done later
  Double_t pard0rphimisal[3]={0,0,0};
  Double_t pard0zMC[3]={85.,130.,0.7};// d0(pt)=[0]+[1]/(pt^[2]); in micron, conversion to cm is done later
  Double_t pard0zmisal[3]={0,0,0};
  if(misal=="data") {
    //use this to reproduce data d0(pt) trend for pions
    pard0rphimisal[0]=37.;
    pard0rphimisal[1]=37.5;
    pard0rphimisal[2]=1.25;
    pard0zmisal[0]=96.;
    pard0zmisal[1]=131.;
    pard0zmisal[2]=0.7;
  }
  else if(misal=="resB") {
    // temporary values: asymptotic value larger by a factor 1.2 w.r.t. MC
    pard0rphimisal[0]=44.4;
    pard0rphimisal[1]=37.5;
    pard0rphimisal[2]=1.25;
    pard0zmisal[0]=115.2;
    pard0zmisal[1]=131.;
    pard0zmisal[2]=0.7;
  }
  else if(misal=="resC") {
    // temporary values: slightly larger asymptotic value, larger values at low pt
    pard0rphimisal[0]=40.;
    pard0rphimisal[1]=40.;
    pard0rphimisal[2]=1.3;
    pard0zmisal[0]=125.;
    pard0zmisal[1]=131.;
    pard0zmisal[2]=0.85;
  }
  else printf("AliAODRecoDecayHF::Misalign():  wrong misalign type specified \n");
 

  AliAODVertex *evVtx=0x0,*secVtx=0x0;
  Double_t evVtxPos[3]={-9999.,-9999.,-9999.},secVtxPos[3]={-9999.,9999.,9999.};
  if(fOwnPrimaryVtx)fOwnPrimaryVtx->GetXYZ(evVtxPos);
  else {
    evVtx=(AliAODVertex*)(fEventPrimaryVtx.GetObject());
    evVtx->GetXYZ(evVtxPos);
  }
  secVtx=(AliAODVertex*)GetSecondaryVtx();
  secVtx->GetXYZ(secVtxPos);
  
  TVector3 v2v1(secVtxPos[0]-evVtxPos[0],secVtxPos[1]-evVtxPos[1],0.);

  Double_t sigmarphinull,sigmarphimisal,sigmarphiadd;
  Double_t sigmaznull,sigmazmisal,sigmazadd;
  Double_t deltad0rphi[10],deltad0z[10];
  
  // loop on the two prongs
  for(Int_t i=0; i<fNProngs; i++) { 
    sigmarphinull = pard0rphiMC[0]+pard0rphiMC[1]/TMath::Power(PtProng(i),pard0rphiMC[2]);
    sigmarphimisal = pard0rphimisal[0]+pard0rphimisal[1]/TMath::Power(PtProng(i),pard0rphimisal[2]);
    if(sigmarphimisal>sigmarphinull) {
      sigmarphiadd = TMath::Sqrt(sigmarphimisal*sigmarphimisal-
				 sigmarphinull*sigmarphinull);
      deltad0rphi[i] = gRandom->Gaus(0.,sigmarphiadd);
    } else {
      deltad0rphi[i] = 0.;
    }

    sigmaznull =  pard0zMC[0]+pard0zMC[1]/TMath::Power(PtProng(i),pard0zMC[2]);
    sigmazmisal = pard0zmisal[0]+pard0zmisal[1]/TMath::Power(PtProng(i),pard0zmisal[2]);
    if(sigmazmisal>sigmaznull) {
      sigmazadd = TMath::Sqrt(sigmazmisal*sigmazmisal-
			      sigmaznull*sigmaznull);
      deltad0z[i] = gRandom->Gaus(0.,sigmazadd);
    } else {
      deltad0z[i] = 0.;
    }

    TVector3 pxy(fPx[i],fPy[i],0.);
    TVector3 pxycrossv2v1=pxy.Cross(v2v1);
    if( pxycrossv2v1.Z()*fd0[i] > 0 ) {
      secVtxPos[0]+=1.e-4*deltad0rphi[i]*(-fPy[i])/PtProng(i);// e-4: conversion to cm
      secVtxPos[1]+=1.e-4*deltad0rphi[i]*(+fPx[i])/PtProng(i);    
    } else {
      secVtxPos[0]+=1.e-4*deltad0rphi[i]*(+fPy[i])/PtProng(i);
      secVtxPos[1]+=1.e-4*deltad0rphi[i]*(-fPx[i])/PtProng(i);    
    }
    
    // change d0rphi
    fd0[i] += 1.e-4*deltad0rphi[i]; // e-4: conversion to cm
    // change secondary vertex z
    secVtxPos[2]+=0.5e-4*deltad0z[i];
  }
  secVtx->SetX(secVtxPos[0]);
  secVtx->SetY(secVtxPos[1]);
  secVtx->SetZ(secVtxPos[2]);

  return;
}
Commit	Line	Data
3244eeed	1	/**************************************************************************
	2	* Copyright(c) 1998-2006, 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
27de2dfb	16	/* $Id$ */
27de2dfb	17
3244eeed	18	/////////////////////////////////////////////////////////////
	19	//
	20	// Base class for AOD reconstructed heavy-flavour decay
	21	//
	22	// Author: A.Dainese, andrea.dainese@lnl.infn.it
	23	/////////////////////////////////////////////////////////////
	24
	25	#include <TDatabasePDG.h>
	26	#include <TVector3.h>
15c6df42	27	#include <TRandom.h>
3244eeed	28	#include "AliAODRecoDecay.h"
3244eeed	29	#include "AliAODRecoDecayHF.h"
ead561f0	30	#include "AliAODEvent.h"
	31	#include "AliVertexerTracks.h"
	32	#include "AliExternalTrackParam.h"
ec653946	33	#include "AliKFVertex.h"
	34	#include "AliVVertex.h"
	35	#include "AliESDVertex.h"
3244eeed	36
	37	ClassImp(AliAODRecoDecayHF)
	38
	39	//--------------------------------------------------------------------------
	40	AliAODRecoDecayHF::AliAODRecoDecayHF() :
	41	AliAODRecoDecay(),
	42	fOwnPrimaryVtx(0x0),
460cd990	43	fEventPrimaryVtx(),
a9b75906	44	fListOfCuts(),
6185d025	45	fd0err(0x0),
939850df	46	fProngID(0x0),
939850df	47	fSelectionMap(0)
3244eeed	48	{
	49	//
	50	// Default Constructor
	51	//
	52	}
	53	//--------------------------------------------------------------------------
	54	AliAODRecoDecayHF::AliAODRecoDecayHF(AliAODVertex *vtx2,Int_t nprongs,Short_t charge,
	55	Double_t px,Double_t py,Double_t *pz,
	56	Double_t d0,Double_t d0err) :
	57	AliAODRecoDecay(vtx2,nprongs,charge,px,py,pz,d0),
	58	fOwnPrimaryVtx(0x0),
460cd990	59	fEventPrimaryVtx(),
a9b75906	60	fListOfCuts(),
6185d025	61	fd0err(0x0),
939850df	62	fProngID(0x0),
939850df	63	fSelectionMap(0)
3244eeed	64	{
	65	//
	66	// Constructor with AliAODVertex for decay vertex
	67	//
	68	fd0err = new Double_t[GetNProngs()];
	69	for(Int_t i=0; i<GetNProngs(); i++) fd0err[i] = d0err[i];
	70	}
	71	//--------------------------------------------------------------------------
	72	AliAODRecoDecayHF::AliAODRecoDecayHF(AliAODVertex *vtx2,Int_t nprongs,Short_t charge,
	73	Double_t d0,Double_t d0err) :
	74	AliAODRecoDecay(vtx2,nprongs,charge,d0),
	75	fOwnPrimaryVtx(0x0),
460cd990	76	fEventPrimaryVtx(),
a9b75906	77	fListOfCuts(),
6185d025	78	fd0err(0x0),
939850df	79	fProngID(0x0),
939850df	80	fSelectionMap(0)
3244eeed	81	{
	82	//
	83	// Constructor with AliAODVertex for decay vertex and without prongs momenta
	84	//
	85	fd0err = new Double_t[GetNProngs()];
	86	for(Int_t i=0; i<GetNProngs(); i++) fd0err[i] = d0err[i];
	87	}
	88	//--------------------------------------------------------------------------
b39168f9	89	AliAODRecoDecayHF::AliAODRecoDecayHF(Double_t vtx1[3],Double_t vtx2[3],
	90	Int_t nprongs,Short_t charge,
	91	Double_t px,Double_t py,Double_t *pz,
	92	Double_t *d0) :
	93	AliAODRecoDecay(0x0,nprongs,charge,px,py,pz,d0),
	94	fOwnPrimaryVtx(0x0),
460cd990	95	fEventPrimaryVtx(),
a9b75906	96	fListOfCuts(),
b39168f9	97	fd0err(0x0),
939850df	98	fProngID(0x0),
939850df	99	fSelectionMap(0)
b39168f9	100	{
	101	//
	102	// Constructor that can used for a "MC" object
	103	//
	104
	105	fOwnPrimaryVtx = new AliAODVertex(vtx1);
	106
	107	AliAODVertex *vtx = new AliAODVertex(vtx2);
	108	SetOwnSecondaryVtx(vtx);
	109
	110	}
	111	//--------------------------------------------------------------------------
3244eeed	112	AliAODRecoDecayHF::AliAODRecoDecayHF(const AliAODRecoDecayHF &source) :
3244eeed	113	AliAODRecoDecay(source),
0a65d33f	114	fOwnPrimaryVtx(0x0),
460cd990	115	fEventPrimaryVtx(source.fEventPrimaryVtx),
a9b75906	116	fListOfCuts(source.fListOfCuts),
6185d025	117	fd0err(0x0),
939850df	118	fProngID(0x0),
939850df	119	fSelectionMap(source.fSelectionMap)
3244eeed	120	{
	121	//
	122	// Copy constructor
	123	//
0a65d33f	124	if(source.GetOwnPrimaryVtx()) fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
0a65d33f	125
3244eeed	126	if(source.GetNProngs()>0) {
	127	fd0err = new Double_t[GetNProngs()];
	128	memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
6185d025	129	if(source.fProngID) {
	130	fProngID = new UShort_t[GetNProngs()];
	131	memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
	132	}
3244eeed	133	}
	134	}
	135	//--------------------------------------------------------------------------
	136	AliAODRecoDecayHF &AliAODRecoDecayHF::operator=(const AliAODRecoDecayHF &source)
	137	{
	138	//
	139	// assignment operator
	140	//
	141	if(&source == this) return *this;
dcb444c9	142
	143	AliAODRecoDecay::operator=(source);
	144
460cd990	145	fEventPrimaryVtx = source.fEventPrimaryVtx;
a9b75906	146	fListOfCuts = source.fListOfCuts;
939850df	147	fSelectionMap = source.fSelectionMap;
460cd990	148
59fee8a6	149	if(source.GetOwnPrimaryVtx()) {
	150	delete fOwnPrimaryVtx;
	151	fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
	152	}
0a65d33f	153
3244eeed	154	if(source.GetNProngs()>0) {
d95be5e9	155	if(source.fd0err) {
59fee8a6	156	delete [] fd0err;
d95be5e9	157	fd0err = new Double_t[GetNProngs()];
	158	memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
	159	}
6185d025	160	if(source.fProngID) {
59fee8a6	161	delete [] fProngID;
6185d025	162	fProngID = new UShort_t[GetNProngs()];
	163	memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
	164	}
3244eeed	165	}
	166	return *this;
	167	}
	168	//--------------------------------------------------------------------------
	169	AliAODRecoDecayHF::~AliAODRecoDecayHF() {
	170	//
	171	// Default Destructor
	172	//
	173	if(fOwnPrimaryVtx) delete fOwnPrimaryVtx;
	174	if(fd0err) delete [] fd0err;
6185d025	175	if(fProngID) delete [] fProngID;
3244eeed	176	}
3244eeed	177	//---------------------------------------------------------------------------
ec653946	178	AliKFParticle AliAODRecoDecayHF::ApplyVertexingKF(Int_t iprongs,Int_t nprongs,Int_t pdgs,Bool_t topoCostraint, Double_t bzkG, Double_t mass) const {
	179	//
	180	// Applies the KF vertexer
	181	// Int_t iprongs[nprongs] = indices of the prongs to be used from the vertexer
	182	// Int_t pdgs[nprongs] = pdgs assigned to the prongs, needed to define the AliKFParticle
	183	// Bool_t topoCostraint = if kTRUE, the topological constraint is applied
	184	// Double_t bzkG = magnetic field
	185	// Double_t mass[2] = {mass, sigma} for the mass constraint (if mass[0]>0 the constraint is applied).
	186	//
	187
	188	AliKFParticle::SetField(bzkG);
	189	AliKFParticle *vertexKF=0;
	190
	191	AliKFVertex copyKF;
	192	Int_t nt=0,ntcheck=0;
	193
	194	Double_t pos[3]={0.,0.,0.};
	195	if(!fOwnPrimaryVtx) {
	196	printf("AliAODRecoDecayHF::ApplyVertexingKF(): cannot apply because primary vertex is not found\n");
	197	return vertexKF;
	198	}
	199	fOwnPrimaryVtx->GetXYZ(pos);
	200	Int_t contr=fOwnPrimaryVtx->GetNContributors();
	201	Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
	202	fOwnPrimaryVtx->GetCovarianceMatrix(covmatrix);
	203	Double_t chi2=fOwnPrimaryVtx->GetChi2();
	204	AliESDVertex primaryVtx2(pos,covmatrix,chi2,contr,"Vertex");
	205
	206
	207	if(topoCostraint){
	208	copyKF=AliKFVertex(primaryVtx2);
	209	nt=primaryVtx2.GetNContributors();
	210	ntcheck=nt;
	211	}
	212
	213	vertexKF = new AliKFParticle();
	214	for(Int_t i= 0;i<nprongs;i++){
	215	Int_t ipr=iprongs[i];
	216	AliAODTrack aodTrack = (AliAODTrack)GetDaughter(ipr);
	217	if(!aodTrack) {
	218	printf("AliAODRecoDecayHF::ApplyVertexingKF(): no daughters available\n");
	219	delete vertexKF; vertexKF=NULL;
	220	return vertexKF;
	221	}
	222	AliKFParticle daughterKF(*aodTrack,pdgs[i]);
	223	vertexKF->AddDaughter(daughterKF);
	224
	225	if(topoCostraint && nt>0){
	226	//Int_t index=(Int_t)GetProngID(ipr);
	227	if(!aodTrack->GetUsedForPrimVtxFit()) continue;
	228	copyKF -= daughterKF;
	229	ntcheck--;
	230	}
	231	}
	232
	233	if(topoCostraint){
	234	if(ntcheck>0) {
	235	copyKF += (*vertexKF);
	236	vertexKF->SetProductionVertex(copyKF);
	237	}
	238	}
	239
	240	if(mass[0]>0.){
	241	vertexKF->SetMassConstraint(mass[0],mass[1]);
242	}
243
244	return vertexKF;
245	}
ead561f0	246	//---------------------------------------------------------------------------
	247	AliAODVertex* AliAODRecoDecayHF::RemoveDaughtersFromPrimaryVtx(AliAODEvent *aod) {
	248	//
	249	// This method returns a primary vertex without the daughter tracks of the
	250	// candidate and it recalculates the impact parameters and errors.
	251	//
	252	// The output vertex is created with "new". The user has to
	253	// set it to the candidate with SetOwnPrimaryVtx(), unset it at the end
	254	// of processing with UnsetOwnPrimaryVtx() and delete it.
	255	// If a NULL pointer is returned, the removal failed (too few tracks left).
	256	//
	257	// For the moment, the primary vertex is recalculated from scratch without
	258	// the daughter tracks.
	259	//
	260
	261	AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
	262	if(!vtxAOD) return 0;
	263	TString title=vtxAOD->GetTitle();
	264	if(!title.Contains("VertexerTracks")) return 0;
	265
	266
	267
	268	AliVertexerTracks *vertexer = new AliVertexerTracks(aod->GetMagneticField());
	269
	270	Int_t ndg = GetNDaughters();
	271
	272	vertexer->SetITSMode();
4a037efc	273	vertexer->SetMinClusters(3);
ead561f0	274	vertexer->SetConstraintOff();
	275
	276	if(title.Contains("WithConstraint")) {
	277	Float_t diamondcovxy[3];
	278	aod->GetDiamondCovXY(diamondcovxy);
	279	Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
	280	Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
	281	AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
	282	vertexer->SetVtxStart(diamond);
	283	delete diamond; diamond=NULL;
	284	}
	285
11d95f54	286	Int_t skipped[10]; for(Int_t i=0;i<10;i++) skipped[i]=-1;
ead561f0	287	Int_t nTrksToSkip=0,id;
	288	AliAODTrack *t = 0;
	289	for(Int_t i=0; i<ndg; i++) {
	290	t = (AliAODTrack*)GetDaughter(i);
	291	id = (Int_t)t->GetID();
	292	if(id<0) continue;
	293	skipped[nTrksToSkip++] = id;
	294	}
	295	vertexer->SetSkipTracks(nTrksToSkip,skipped);
	296	AliESDVertex *vtxESDNew = vertexer->FindPrimaryVertex(aod);
	297
	298	delete vertexer; vertexer=NULL;
	299
	300	if(!vtxESDNew) return 0;
	301	if(vtxESDNew->GetNContributors()<=0) {
	302	delete vtxESDNew; vtxESDNew=NULL;
	303	return 0;
	304	}
	305
	306	// convert to AliAODVertex
	307	Double_t pos[3],cov[6],chi2perNDF;
	308	vtxESDNew->GetXYZ(pos); // position
	309	vtxESDNew->GetCovMatrix(cov); //covariance matrix
	310	chi2perNDF = vtxESDNew->GetChi2toNDF();
	311	delete vtxESDNew; vtxESDNew=NULL;
	312
	313	AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);
	314
70841bf2	315	RecalculateImpPars(vtxAODNew,aod);
	316
	317	return vtxAODNew;
	318	}
	319	//-----------------------------------------------------------------------------------
	320	void AliAODRecoDecayHF::RecalculateImpPars(AliAODVertex vtxAODNew,AliAODEvent aod) {
	321	//
ead561f0	322	// now recalculate the daughters impact parameters
70841bf2	323	//
ead561f0	324	Double_t dz[2],covdz[3];
70841bf2	325	for(Int_t i=0; i<GetNDaughters(); i++) {
70841bf2	326	AliAODTrack t = (AliAODTrack)GetDaughter(i);
fb690bb8	327	AliExternalTrackParam etp; etp.CopyFromVTrack(t);
4ec1ca0f	328	if(etp.PropagateToDCA(vtxAODNew,aod->GetMagneticField(),3.,dz,covdz)) {
ead561f0	329	fd0[i] = dz[0];
	330	fd0err[i] = TMath::Sqrt(covdz[0]);
	331	}
ead561f0	332	}
ead561f0	333
70841bf2	334	return;
ead561f0	335	}
15c6df42	336	//-----------------------------------------------------------------------------------
	337	void AliAODRecoDecayHF::Misalign(TString misal) {
	338	//
	339	// Method to smear the impact parameter of the duaghter tracks
	340	// and the sec. vtx position accordinlgy
	341	// Useful to study effect of misalignment.
	342	// The starting point are parameterisations of the impact parameter resolution
	343	// from MC and data
	344	// Errors on d0 and vtx are not recalculated (to be done)
	345	//
	346	if(misal=="null")return;
	347	Double_t pard0rphiMC[3]={36.7,36.,1.25};// d0(pt)=[0]+[1]/(pt^[2]); in micron, conversion to cm is done later
	348	Double_t pard0rphimisal[3]={0,0,0};
	349	Double_t pard0zMC[3]={85.,130.,0.7};// d0(pt)=[0]+[1]/(pt^[2]); in micron, conversion to cm is done later
	350	Double_t pard0zmisal[3]={0,0,0};
	351	if(misal=="data") {
	352	//use this to reproduce data d0(pt) trend for pions
	353	pard0rphimisal[0]=37.;
	354	pard0rphimisal[1]=37.5;
	355	pard0rphimisal[2]=1.25;
	356	pard0zmisal[0]=96.;
	357	pard0zmisal[1]=131.;
	358	pard0zmisal[2]=0.7;
	359	}
	360	else if(misal=="resB") {
	361	// temporary values: asymptotic value larger by a factor 1.2 w.r.t. MC
	362	pard0rphimisal[0]=44.4;
	363	pard0rphimisal[1]=37.5;
	364	pard0rphimisal[2]=1.25;
	365	pard0zmisal[0]=115.2;
	366	pard0zmisal[1]=131.;
	367	pard0zmisal[2]=0.7;
	368	}
	369	else if(misal=="resC") {
	370	// temporary values: slightly larger asymptotic value, larger values at low pt
	371	pard0rphimisal[0]=40.;
	372	pard0rphimisal[1]=40.;
	373	pard0rphimisal[2]=1.3;
	374	pard0zmisal[0]=125.;
	375	pard0zmisal[1]=131.;
	376	pard0zmisal[2]=0.85;
	377	}
	378	else printf("AliAODRecoDecayHF::Misalign(): wrong misalign type specified \n");
	379
	380
	381	AliAODVertex evVtx=0x0,secVtx=0x0;
	382	Double_t evVtxPos[3]={-9999.,-9999.,-9999.},secVtxPos[3]={-9999.,9999.,9999.};
	383	if(fOwnPrimaryVtx)fOwnPrimaryVtx->GetXYZ(evVtxPos);
	384	else {
	385	evVtx=(AliAODVertex*)(fEventPrimaryVtx.GetObject());
	386	evVtx->GetXYZ(evVtxPos);
	387	}
	388	secVtx=(AliAODVertex*)GetSecondaryVtx();
	389	secVtx->GetXYZ(secVtxPos);
	390
	391	TVector3 v2v1(secVtxPos[0]-evVtxPos[0],secVtxPos[1]-evVtxPos[1],0.);
	392
	393	Double_t sigmarphinull,sigmarphimisal,sigmarphiadd;
	394	Double_t sigmaznull,sigmazmisal,sigmazadd;
	395	Double_t deltad0rphi[10],deltad0z[10];
	396
	397	// loop on the two prongs
	398	for(Int_t i=0; i<fNProngs; i++) {
	399	sigmarphinull = pard0rphiMC[0]+pard0rphiMC[1]/TMath::Power(PtProng(i),pard0rphiMC[2]);
400	sigmarphimisal = pard0rphimisal[0]+pard0rphimisal[1]/TMath::Power(PtProng(i),pard0rphimisal[2]);
401	if(sigmarphimisal>sigmarphinull) {
402	sigmarphiadd = TMath::Sqrt(sigmarphimisal*sigmarphimisal-
403	sigmarphinull*sigmarphinull);
404	deltad0rphi[i] = gRandom->Gaus(0.,sigmarphiadd);
405	} else {
406	deltad0rphi[i] = 0.;
407	}
408
409	sigmaznull = pard0zMC[0]+pard0zMC[1]/TMath::Power(PtProng(i),pard0zMC[2]);
410	sigmazmisal = pard0zmisal[0]+pard0zmisal[1]/TMath::Power(PtProng(i),pard0zmisal[2]);
411	if(sigmazmisal>sigmaznull) {
412	sigmazadd = TMath::Sqrt(sigmazmisal*sigmazmisal-
413	sigmaznull*sigmaznull);
414	deltad0z[i] = gRandom->Gaus(0.,sigmazadd);
415	} else {
416	deltad0z[i] = 0.;
417	}
418
419	TVector3 pxy(fPx[i],fPy[i],0.);
420	TVector3 pxycrossv2v1=pxy.Cross(v2v1);
421	if( pxycrossv2v1.Z()*fd0[i] > 0 ) {
422	secVtxPos[0]+=1.e-4deltad0rphi[i](-fPy[i])/PtProng(i);// e-4: conversion to cm
423	secVtxPos[1]+=1.e-4deltad0rphi[i](+fPx[i])/PtProng(i);
424	} else {
425	secVtxPos[0]+=1.e-4deltad0rphi[i](+fPy[i])/PtProng(i);
426	secVtxPos[1]+=1.e-4deltad0rphi[i](-fPx[i])/PtProng(i);
427	}
428
429	// change d0rphi
430	fd0[i] += 1.e-4*deltad0rphi[i]; // e-4: conversion to cm
431	// change secondary vertex z
432	secVtxPos[2]+=0.5e-4*deltad0z[i];
433	}
434	secVtx->SetX(secVtxPos[0]);
435	secVtx->SetY(secVtxPos[1]);
436	secVtx->SetZ(secVtxPos[2]);
437
438	return;
439	}