Added LHC10h run list for flow analysis (Giacomo)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAODRecoDecayHF.cxx
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
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>
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
37 ClassImp(AliAODRecoDecayHF)
38
39 //--------------------------------------------------------------------------
40 AliAODRecoDecayHF::AliAODRecoDecayHF() :
41   AliAODRecoDecay(),
42   fOwnPrimaryVtx(0x0),
43   fEventPrimaryVtx(),
44   fListOfCuts(),
45   fd0err(0x0), 
46   fProngID(0x0),
47   fSelectionMap(0)
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),
59   fEventPrimaryVtx(),
60   fListOfCuts(),
61   fd0err(0x0),
62   fProngID(0x0),
63   fSelectionMap(0)
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),
76   fEventPrimaryVtx(),
77   fListOfCuts(),
78   fd0err(0x0),
79   fProngID(0x0),
80   fSelectionMap(0)
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 //--------------------------------------------------------------------------
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),
95   fEventPrimaryVtx(),
96   fListOfCuts(),
97   fd0err(0x0),
98   fProngID(0x0), 
99   fSelectionMap(0)
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 //--------------------------------------------------------------------------
112 AliAODRecoDecayHF::AliAODRecoDecayHF(const AliAODRecoDecayHF &source) :
113   AliAODRecoDecay(source),
114   fOwnPrimaryVtx(0x0),
115   fEventPrimaryVtx(source.fEventPrimaryVtx),
116   fListOfCuts(source.fListOfCuts),
117   fd0err(0x0),
118   fProngID(0x0),
119   fSelectionMap(source.fSelectionMap)
120 {
121   //
122   // Copy constructor
123   //
124   if(source.GetOwnPrimaryVtx()) fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
125
126   if(source.GetNProngs()>0) {
127     fd0err = new Double_t[GetNProngs()];
128     memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
129     if(source.fProngID) {
130       fProngID = new UShort_t[GetNProngs()];
131       memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
132     }
133   }
134 }
135 //--------------------------------------------------------------------------
136 AliAODRecoDecayHF &AliAODRecoDecayHF::operator=(const AliAODRecoDecayHF &source)
137 {
138   //
139   // assignment operator
140   //
141   if(&source == this) return *this;
142
143   AliAODRecoDecay::operator=(source);
144
145   fEventPrimaryVtx = source.fEventPrimaryVtx;
146   fListOfCuts = source.fListOfCuts;
147   fSelectionMap = source.fSelectionMap;
148
149   if(source.GetOwnPrimaryVtx()) fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
150
151   if(source.GetNProngs()>0) {
152     fd0err = new Double_t[GetNProngs()];
153     memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
154     if(source.fProngID) {
155       fProngID = new UShort_t[GetNProngs()];
156       memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
157     }
158   }
159   return *this;
160 }
161 //--------------------------------------------------------------------------
162 AliAODRecoDecayHF::~AliAODRecoDecayHF() {
163   //  
164   // Default Destructor
165   //
166   if(fOwnPrimaryVtx) delete fOwnPrimaryVtx;
167   if(fd0err) delete [] fd0err;
168   if(fProngID) delete [] fProngID;
169 }
170 //---------------------------------------------------------------------------
171 AliKFParticle *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 }
239 //---------------------------------------------------------------------------
240 AliAODVertex* 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();
266   vertexer->SetMinClusters(3);
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
308   RecalculateImpPars(vtxAODNew,aod);
309
310   return vtxAODNew;
311 }
312 //-----------------------------------------------------------------------------------
313 void AliAODRecoDecayHF::RecalculateImpPars(AliAODVertex *vtxAODNew,AliAODEvent* aod) {
314   //
315   // now recalculate the daughters impact parameters
316   //
317   AliExternalTrackParam *etp = 0;
318   Double_t dz[2],covdz[3];
319   for(Int_t i=0; i<GetNDaughters(); i++) {
320     AliAODTrack *t = (AliAODTrack*)GetDaughter(i);
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
329   return;
330 }
331 //-----------------------------------------------------------------------------------
332 void 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 }