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