Added method to recalculate the primary vertex without the daughter tracks
[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 /////////////////////////////////////////////////////////////
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"
27 #include "AliAODEvent.h"
28 #include "AliVertexerTracks.h"
29 #include "AliExternalTrackParam.h"
30 #include "AliKFVertex.h"
31 #include "AliVVertex.h"
32 #include "AliESDVertex.h"
33
34 ClassImp(AliAODRecoDecayHF)
35
36 //--------------------------------------------------------------------------
37 AliAODRecoDecayHF::AliAODRecoDecayHF() :
38   AliAODRecoDecay(),
39   fOwnPrimaryVtx(0x0),
40   fEventPrimaryVtx(),
41   fListOfCuts(),
42   fd0err(0x0), 
43   fProngID(0x0) 
44 {
45   //
46   // Default Constructor
47   //
48 }
49 //--------------------------------------------------------------------------
50 AliAODRecoDecayHF::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),
55   fEventPrimaryVtx(),
56   fListOfCuts(),
57   fd0err(0x0),
58   fProngID(0x0) 
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 //--------------------------------------------------------------------------
67 AliAODRecoDecayHF::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),
71   fEventPrimaryVtx(),
72   fListOfCuts(),
73   fd0err(0x0),
74   fProngID(0x0) 
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 //--------------------------------------------------------------------------
83 AliAODRecoDecayHF::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),
89   fEventPrimaryVtx(),
90   fListOfCuts(),
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 //--------------------------------------------------------------------------
105 AliAODRecoDecayHF::AliAODRecoDecayHF(const AliAODRecoDecayHF &source) :
106   AliAODRecoDecay(source),
107   fOwnPrimaryVtx(0x0),
108   fEventPrimaryVtx(source.fEventPrimaryVtx),
109   fListOfCuts(source.fListOfCuts),
110   fd0err(0x0),
111   fProngID(0x0)
112 {
113   //
114   // Copy constructor
115   //
116   if(source.GetOwnPrimaryVtx()) fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
117
118   if(source.GetNProngs()>0) {
119     fd0err = new Double_t[GetNProngs()];
120     memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
121     if(source.fProngID) {
122       fProngID = new UShort_t[GetNProngs()];
123       memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
124     }
125   }
126 }
127 //--------------------------------------------------------------------------
128 AliAODRecoDecayHF &AliAODRecoDecayHF::operator=(const AliAODRecoDecayHF &source)
129 {
130   //
131   // assignment operator
132   //
133   if(&source == this) return *this;
134
135   AliAODRecoDecay::operator=(source);
136
137   fEventPrimaryVtx = source.fEventPrimaryVtx;
138   fListOfCuts = source.fListOfCuts;
139
140   if(source.GetOwnPrimaryVtx()) fOwnPrimaryVtx = new AliAODVertex(*(source.GetOwnPrimaryVtx()));
141
142   if(source.GetNProngs()>0) {
143     fd0err = new Double_t[GetNProngs()];
144     memcpy(fd0err,source.fd0err,GetNProngs()*sizeof(Double_t));
145     if(source.fProngID) {
146       fProngID = new UShort_t[GetNProngs()];
147       memcpy(fProngID,source.fProngID,GetNProngs()*sizeof(UShort_t));
148     }
149   }
150   return *this;
151 }
152 //--------------------------------------------------------------------------
153 AliAODRecoDecayHF::~AliAODRecoDecayHF() {
154   //  
155   // Default Destructor
156   //
157   if(fOwnPrimaryVtx) delete fOwnPrimaryVtx;
158   if(fd0err) delete [] fd0err;
159   if(fProngID) delete [] fProngID;
160 }
161 //---------------------------------------------------------------------------
162 AliKFParticle *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 }
230 //---------------------------------------------------------------------------
231 AliAODVertex* 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 }