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