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