Added method to recalculate the primary vertex without the daughter tracks
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAODRecoDecayHF.cxx
CommitLineData
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
16/////////////////////////////////////////////////////////////
17//
18// Base class for AOD reconstructed heavy-flavour decay
19//
20// Author: A.Dainese, andrea.dainese@lnl.infn.it
21/////////////////////////////////////////////////////////////
22
23#include <TDatabasePDG.h>
24#include <TVector3.h>
25#include "AliAODRecoDecay.h"
26#include "AliAODRecoDecayHF.h"
ead561f0 27#include "AliAODEvent.h"
28#include "AliVertexerTracks.h"
29#include "AliExternalTrackParam.h"
ec653946 30#include "AliKFVertex.h"
31#include "AliVVertex.h"
32#include "AliESDVertex.h"
3244eeed 33
34ClassImp(AliAODRecoDecayHF)
35
36//--------------------------------------------------------------------------
37AliAODRecoDecayHF::AliAODRecoDecayHF() :
38 AliAODRecoDecay(),
39 fOwnPrimaryVtx(0x0),
460cd990 40 fEventPrimaryVtx(),
a9b75906 41 fListOfCuts(),
6185d025 42 fd0err(0x0),
43 fProngID(0x0)
3244eeed 44{
45 //
46 // Default Constructor
47 //
48}
49//--------------------------------------------------------------------------
50AliAODRecoDecayHF::AliAODRecoDecayHF(AliAODVertex *vtx2,Int_t nprongs,Short_t charge,
51 Double_t *px,Double_t *py,Double_t *pz,
52 Double_t *d0,Double_t *d0err) :
53 AliAODRecoDecay(vtx2,nprongs,charge,px,py,pz,d0),
54 fOwnPrimaryVtx(0x0),
460cd990 55 fEventPrimaryVtx(),
a9b75906 56 fListOfCuts(),
6185d025 57 fd0err(0x0),
58 fProngID(0x0)
3244eeed 59{
60 //
61 // Constructor with AliAODVertex for decay vertex
62 //
63 fd0err = new Double_t[GetNProngs()];
64 for(Int_t i=0; i<GetNProngs(); i++) fd0err[i] = d0err[i];
65}
66//--------------------------------------------------------------------------
67AliAODRecoDecayHF::AliAODRecoDecayHF(AliAODVertex *vtx2,Int_t nprongs,Short_t charge,
68 Double_t *d0,Double_t *d0err) :
69 AliAODRecoDecay(vtx2,nprongs,charge,d0),
70 fOwnPrimaryVtx(0x0),
460cd990 71 fEventPrimaryVtx(),
a9b75906 72 fListOfCuts(),
6185d025 73 fd0err(0x0),
74 fProngID(0x0)
3244eeed 75{
76 //
77 // Constructor with AliAODVertex for decay vertex and without prongs momenta
78 //
79 fd0err = new Double_t[GetNProngs()];
80 for(Int_t i=0; i<GetNProngs(); i++) fd0err[i] = d0err[i];
81}
82//--------------------------------------------------------------------------
b39168f9 83AliAODRecoDecayHF::AliAODRecoDecayHF(Double_t vtx1[3],Double_t vtx2[3],
84 Int_t nprongs,Short_t charge,
85 Double_t *px,Double_t *py,Double_t *pz,
86 Double_t *d0) :
87 AliAODRecoDecay(0x0,nprongs,charge,px,py,pz,d0),
88 fOwnPrimaryVtx(0x0),
460cd990 89 fEventPrimaryVtx(),
a9b75906 90 fListOfCuts(),
b39168f9 91 fd0err(0x0),
92 fProngID(0x0)
93{
94 //
95 // Constructor that can used for a "MC" object
96 //
97
98 fOwnPrimaryVtx = new AliAODVertex(vtx1);
99
100 AliAODVertex *vtx = new AliAODVertex(vtx2);
101 SetOwnSecondaryVtx(vtx);
102
103}
104//--------------------------------------------------------------------------
3244eeed 105AliAODRecoDecayHF::AliAODRecoDecayHF(const AliAODRecoDecayHF &source) :
106 AliAODRecoDecay(source),
0a65d33f 107 fOwnPrimaryVtx(0x0),
460cd990 108 fEventPrimaryVtx(source.fEventPrimaryVtx),
a9b75906 109 fListOfCuts(source.fListOfCuts),
6185d025 110 fd0err(0x0),
111 fProngID(0x0)
3244eeed 112{
113 //
114 // Copy constructor
115 //
0a65d33f 116 if(source.GetOwnPrimaryVtx()) fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
117
3244eeed 118 if(source.GetNProngs()>0) {
119 fd0err = new Double_t[GetNProngs()];
120 memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
6185d025 121 if(source.fProngID) {
122 fProngID = new UShort_t[GetNProngs()];
123 memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
124 }
3244eeed 125 }
126}
127//--------------------------------------------------------------------------
128AliAODRecoDecayHF &AliAODRecoDecayHF::operator=(const AliAODRecoDecayHF &source)
129{
130 //
131 // assignment operator
132 //
133 if(&source == this) return *this;
dcb444c9 134
135 AliAODRecoDecay::operator=(source);
136
460cd990 137 fEventPrimaryVtx = source.fEventPrimaryVtx;
a9b75906 138 fListOfCuts = source.fListOfCuts;
460cd990 139
0a65d33f 140 if(source.GetOwnPrimaryVtx()) fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
141
3244eeed 142 if(source.GetNProngs()>0) {
3244eeed 143 fd0err = new Double_t[GetNProngs()];
3244eeed 144 memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
6185d025 145 if(source.fProngID) {
146 fProngID = new UShort_t[GetNProngs()];
147 memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
148 }
3244eeed 149 }
150 return *this;
151}
152//--------------------------------------------------------------------------
153AliAODRecoDecayHF::~AliAODRecoDecayHF() {
154 //
155 // Default Destructor
156 //
157 if(fOwnPrimaryVtx) delete fOwnPrimaryVtx;
158 if(fd0err) delete [] fd0err;
6185d025 159 if(fProngID) delete [] fProngID;
3244eeed 160}
161//---------------------------------------------------------------------------
ec653946 162AliKFParticle *AliAODRecoDecayHF::ApplyVertexingKF(Int_t *iprongs,Int_t nprongs,Int_t *pdgs,Bool_t topoCostraint, Double_t bzkG, Double_t *mass) const {
163 //
164 // Applies the KF vertexer
165 // Int_t iprongs[nprongs] = indices of the prongs to be used from the vertexer
166 // Int_t pdgs[nprongs] = pdgs assigned to the prongs, needed to define the AliKFParticle
167 // Bool_t topoCostraint = if kTRUE, the topological constraint is applied
168 // Double_t bzkG = magnetic field
169 // Double_t mass[2] = {mass, sigma} for the mass constraint (if mass[0]>0 the constraint is applied).
170 //
171
172 AliKFParticle::SetField(bzkG);
173 AliKFParticle *vertexKF=0;
174
175 AliKFVertex copyKF;
176 Int_t nt=0,ntcheck=0;
177
178 Double_t pos[3]={0.,0.,0.};
179 if(!fOwnPrimaryVtx) {
180 printf("AliAODRecoDecayHF::ApplyVertexingKF(): cannot apply because primary vertex is not found\n");
181 return vertexKF;
182 }
183 fOwnPrimaryVtx->GetXYZ(pos);
184 Int_t contr=fOwnPrimaryVtx->GetNContributors();
185 Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
186 fOwnPrimaryVtx->GetCovarianceMatrix(covmatrix);
187 Double_t chi2=fOwnPrimaryVtx->GetChi2();
188 AliESDVertex primaryVtx2(pos,covmatrix,chi2,contr,"Vertex");
189
190
191 if(topoCostraint){
192 copyKF=AliKFVertex(primaryVtx2);
193 nt=primaryVtx2.GetNContributors();
194 ntcheck=nt;
195 }
196
197 vertexKF = new AliKFParticle();
198 for(Int_t i= 0;i<nprongs;i++){
199 Int_t ipr=iprongs[i];
200 AliAODTrack *aodTrack = (AliAODTrack*)GetDaughter(ipr);
201 if(!aodTrack) {
202 printf("AliAODRecoDecayHF::ApplyVertexingKF(): no daughters available\n");
203 delete vertexKF; vertexKF=NULL;
204 return vertexKF;
205 }
206 AliKFParticle daughterKF(*aodTrack,pdgs[i]);
207 vertexKF->AddDaughter(daughterKF);
208
209 if(topoCostraint && nt>0){
210 //Int_t index=(Int_t)GetProngID(ipr);
211 if(!aodTrack->GetUsedForPrimVtxFit()) continue;
212 copyKF -= daughterKF;
213 ntcheck--;
214 }
215 }
216
217 if(topoCostraint){
218 if(ntcheck>0) {
219 copyKF += (*vertexKF);
220 vertexKF->SetProductionVertex(copyKF);
221 }
222 }
223
224 if(mass[0]>0.){
225 vertexKF->SetMassConstraint(mass[0],mass[1]);
226 }
227
228 return vertexKF;
229}
ead561f0 230//---------------------------------------------------------------------------
231AliAODVertex* AliAODRecoDecayHF::RemoveDaughtersFromPrimaryVtx(AliAODEvent *aod) {
232 //
233 // This method returns a primary vertex without the daughter tracks of the
234 // candidate and it recalculates the impact parameters and errors.
235 //
236 // The output vertex is created with "new". The user has to
237 // set it to the candidate with SetOwnPrimaryVtx(), unset it at the end
238 // of processing with UnsetOwnPrimaryVtx() and delete it.
239 // If a NULL pointer is returned, the removal failed (too few tracks left).
240 //
241 // For the moment, the primary vertex is recalculated from scratch without
242 // the daughter tracks.
243 //
244
245 AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
246 if(!vtxAOD) return 0;
247 TString title=vtxAOD->GetTitle();
248 if(!title.Contains("VertexerTracks")) return 0;
249
250
251
252 AliVertexerTracks *vertexer = new AliVertexerTracks(aod->GetMagneticField());
253
254 Int_t ndg = GetNDaughters();
255
256 vertexer->SetITSMode();
257 vertexer->SetMinClusters(4);
258 vertexer->SetConstraintOff();
259
260 if(title.Contains("WithConstraint")) {
261 Float_t diamondcovxy[3];
262 aod->GetDiamondCovXY(diamondcovxy);
263 Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
264 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
265 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
266 vertexer->SetVtxStart(diamond);
267 delete diamond; diamond=NULL;
268 }
269
270 Int_t skipped[10];
271 Int_t nTrksToSkip=0,id;
272 AliAODTrack *t = 0;
273 for(Int_t i=0; i<ndg; i++) {
274 t = (AliAODTrack*)GetDaughter(i);
275 id = (Int_t)t->GetID();
276 if(id<0) continue;
277 skipped[nTrksToSkip++] = id;
278 }
279 vertexer->SetSkipTracks(nTrksToSkip,skipped);
280 AliESDVertex *vtxESDNew = vertexer->FindPrimaryVertex(aod);
281
282 delete vertexer; vertexer=NULL;
283
284 if(!vtxESDNew) return 0;
285 if(vtxESDNew->GetNContributors()<=0) {
286 delete vtxESDNew; vtxESDNew=NULL;
287 return 0;
288 }
289
290 // convert to AliAODVertex
291 Double_t pos[3],cov[6],chi2perNDF;
292 vtxESDNew->GetXYZ(pos); // position
293 vtxESDNew->GetCovMatrix(cov); //covariance matrix
294 chi2perNDF = vtxESDNew->GetChi2toNDF();
295 delete vtxESDNew; vtxESDNew=NULL;
296
297 AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);
298
299 // now recalculate the daughters impact parameters
300 AliExternalTrackParam *etp = 0;
301 Double_t dz[2],covdz[3];
302 for(Int_t i=0; i<ndg; i++) {
303 t = (AliAODTrack*)GetDaughter(i);
304 etp = new AliExternalTrackParam(t);
305 if(etp->PropagateToDCA(vtxAODNew,aod->GetMagneticField(),3.,dz,covdz)) {
306 fd0[i] = dz[0];
307 fd0err[i] = TMath::Sqrt(covdz[0]);
308 }
309 delete etp; etp=NULL;
310 }
311
312 return vtxAODNew;
313}