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