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