]>
Commit | Line | Data |
---|---|---|
a6c5a2e9 | 1 | /************************************************************************** |
2 | * Copyright(c) 2007-2009, 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 | #include <TMath.h> | |
686e4710 | 17 | #include <TRandom.h> |
b176db3b | 18 | #include <TProfile.h> |
19 | #include <TClonesArray.h> | |
20 | #include <TH1F.h> | |
21 | #include <TH2F.h> | |
22 | #include <TF1.h> | |
a6c5a2e9 | 23 | #include "AliAODEvent.h" |
b176db3b | 24 | #include "AliAODMCHeader.h" |
7f51a627 | 25 | #include "AliGenEventHeader.h" |
25504150 | 26 | #include "AliAODMCParticle.h" |
b176db3b | 27 | #include "AliAODRecoDecayHF.h" |
a6c5a2e9 | 28 | #include "AliVertexingHFUtils.h" |
29 | ||
d947ea9c | 30 | /* $Id$ */ |
a6c5a2e9 | 31 | |
32 | /////////////////////////////////////////////////////////////////// | |
33 | // // | |
34 | // Class with functions useful for different D2H analyses // | |
35 | // - event plane resolution // | |
36 | // - <pt> calculation with side band subtraction // | |
25504150 | 37 | // - tracklet multiplicity calculation // |
a6c5a2e9 | 38 | // Origin: F.Prino, Torino, prino@to.infn.it // |
39 | // // | |
40 | /////////////////////////////////////////////////////////////////// | |
41 | ||
42 | ClassImp(AliVertexingHFUtils) | |
43 | ||
44 | //______________________________________________________________________ | |
45 | AliVertexingHFUtils::AliVertexingHFUtils():TObject(), | |
46 | fK(1), | |
47 | fSubRes(1.), | |
48 | fMinEtaForTracklets(-1.), | |
49 | fMaxEtaForTracklets(1.) | |
50 | { | |
51 | // Default contructor | |
52 | } | |
53 | ||
54 | ||
55 | //______________________________________________________________________ | |
56 | AliVertexingHFUtils::AliVertexingHFUtils(Int_t k): | |
57 | TObject(), | |
58 | fK(k), | |
59 | fSubRes(1.), | |
60 | fMinEtaForTracklets(-1.), | |
61 | fMaxEtaForTracklets(1.) | |
62 | { | |
63 | // Standard constructor | |
64 | } | |
65 | ||
66 | ||
5d1396f9 | 67 | //______________________________________________________________________ |
68 | void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance,Double_t &errsignificance){ | |
69 | // calculate significance from S, B and errors | |
70 | ||
71 | ||
72 | Double_t errSigSq=errsignal*errsignal; | |
73 | Double_t errBkgSq=errbackground*errbackground; | |
74 | Double_t sigPlusBkg=signal+background; | |
75 | if(sigPlusBkg>0. && signal>0.){ | |
76 | significance = signal/TMath::Sqrt(signal+background); | |
77 | errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal); | |
78 | }else{ | |
79 | significance=0.; | |
80 | errsignificance=0.; | |
81 | } | |
82 | return; | |
83 | ||
84 | } | |
a6c5a2e9 | 85 | //______________________________________________________________________ |
86 | Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){ | |
87 | // compute chi from polynomial approximation | |
88 | Double_t c[5]; | |
89 | if(k==1){ | |
90 | c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283; | |
91 | } | |
92 | else if(k==2){ | |
93 | c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815; | |
94 | } | |
95 | else return -1; | |
96 | return c[0]*x+c[1]*x*x+c[2]*x*x*x+c[3]*x*x*x*x+c[4]*x*x*x*x*x; | |
97 | } | |
98 | ||
99 | //______________________________________________________________________ | |
100 | Double_t AliVertexingHFUtils:: ResolK1(Double_t x){ | |
101 | return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4)); | |
102 | } | |
103 | ||
104 | ||
105 | //______________________________________________________________________ | |
106 | Double_t AliVertexingHFUtils::FindChi(Double_t res, Int_t k){ | |
107 | // compute chi variable (=v2*sqrt(N)) from external values | |
108 | ||
109 | Double_t x1=0; | |
110 | Double_t x2=15; | |
111 | Double_t y1,y2; | |
112 | if(k==1){ | |
113 | y1=ResolK1(x1)-res; | |
114 | y2=ResolK1(x2)-res; | |
115 | } | |
116 | else if(k==2){ | |
117 | y1=Pol(x1,2)-res; | |
118 | y2=Pol(x2,2)-res; | |
119 | } | |
120 | else return -1; | |
121 | ||
122 | if(y1*y2>0) return -1; | |
123 | if(y1==0) return y1; | |
124 | if(y2==0) return y2; | |
125 | Double_t xmed,ymed; | |
126 | Int_t jiter=0; | |
127 | while((x2-x1)>0.0001){ | |
128 | xmed=0.5*(x1+x2); | |
129 | if(k==1){ | |
130 | y1=ResolK1(x1)-res; | |
131 | y2=ResolK1(x2)-res; | |
132 | ymed=ResolK1(xmed)-res; | |
133 | } | |
134 | else if(k==2){ | |
135 | y1=Pol(x1,2)-res; | |
136 | y2=Pol(x2,2)-res; | |
137 | ymed=Pol(xmed,2)-res; | |
138 | } | |
139 | else return -1; | |
140 | if((y1*ymed)<0) x2=xmed; | |
141 | if((y2*ymed)<0)x1=xmed; | |
142 | if(ymed==0) return xmed; | |
143 | jiter++; | |
144 | } | |
145 | return 0.5*(x1+x2); | |
146 | } | |
147 | ||
148 | //______________________________________________________________________ | |
149 | Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){ | |
150 | // computes event plane resolution starting from sub event resolution | |
151 | Double_t chisub=FindChi(resSub,k); | |
152 | Double_t chifull=chisub*TMath::Sqrt(2); | |
153 | if(k==1) return ResolK1(chifull); | |
154 | else if(k==2) return Pol(chifull,2); | |
155 | else{ | |
156 | printf("k should be <=2\n"); | |
157 | return 1.; | |
158 | } | |
159 | } | |
160 | ||
161 | //______________________________________________________________________ | |
162 | Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){ | |
163 | // computes event plane resolution starting from sub event correlation histogram | |
164 | if(!hSubEvCorr) return 1.; | |
165 | Double_t resSub=GetSubEvResol(hSubEvCorr); | |
166 | return GetFullEvResol(resSub,k); | |
167 | } | |
168 | //______________________________________________________________________ | |
169 | Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){ | |
170 | // computes low limit event plane resolution starting from sub event correlation histogram | |
171 | if(!hSubEvCorr) return 1.; | |
172 | Double_t resSub=GetSubEvResolLowLim(hSubEvCorr); | |
173 | printf("%f\n",resSub); | |
174 | return GetFullEvResol(resSub,k); | |
175 | } | |
176 | //______________________________________________________________________ | |
177 | Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){ | |
178 | // computes high limit event plane resolution starting from sub event correlation histogram | |
179 | if(!hSubEvCorr) return 1.; | |
180 | Double_t resSub=GetSubEvResolHighLim(hSubEvCorr); | |
181 | printf("%f\n",resSub); | |
182 | return GetFullEvResol(resSub,k); | |
183 | } | |
184 | //______________________________________________________________________ | |
185 | Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){ | |
186 | // counts tracklets in given eta range | |
187 | AliAODTracklets* tracklets=ev->GetTracklets(); | |
188 | Int_t nTr=tracklets->GetNumberOfTracklets(); | |
189 | Int_t count=0; | |
190 | for(Int_t iTr=0; iTr<nTr; iTr++){ | |
191 | Double_t theta=tracklets->GetTheta(iTr); | |
192 | Double_t eta=-TMath::Log(TMath::Tan(theta/2.)); | |
193 | if(eta>mineta && eta<maxeta) count++; | |
194 | } | |
195 | return count; | |
196 | } | |
86d84fa4 | 197 | //______________________________________________________________________ |
25504150 | 198 | Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){ |
199 | // counts generated particles in fgiven eta range | |
200 | ||
201 | Int_t nChargedMC=0; | |
202 | for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){ | |
203 | AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i); | |
204 | Int_t charge = part->Charge(); | |
205 | Double_t eta = part->Eta(); | |
206 | if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++; | |
207 | } | |
208 | return nChargedMC; | |
209 | } | |
210 | //______________________________________________________________________ | |
211 | Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){ | |
212 | // counts generated primary particles in given eta range | |
213 | ||
214 | Int_t nChargedMC=0; | |
215 | for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){ | |
216 | AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i); | |
217 | Int_t charge = part->Charge(); | |
218 | Double_t eta = part->Eta(); | |
219 | if(charge!=0 && eta>mineta && eta<maxeta){ | |
220 | if(part->IsPrimary())nChargedMC++; | |
221 | } | |
222 | } | |
223 | return nChargedMC; | |
224 | } | |
225 | //______________________________________________________________________ | |
226 | Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){ | |
227 | // counts generated primary particles in given eta range | |
228 | ||
229 | Int_t nChargedMC=0; | |
230 | for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){ | |
231 | AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i); | |
232 | Int_t charge = part->Charge(); | |
233 | Double_t eta = part->Eta(); | |
234 | if(charge!=0 && eta>mineta && eta<maxeta){ | |
235 | if(part->IsPhysicalPrimary())nChargedMC++; | |
236 | } | |
237 | } | |
238 | return nChargedMC; | |
239 | } | |
240 | //______________________________________________________________________ | |
19df7ca5 | 241 | void AliVertexingHFUtils::AveragePt(Float_t& averagePt, Float_t& errorPt,Float_t ptmin,Float_t ptmax, TH2F* hMassD, Float_t massFromFit, Float_t sigmaFromFit, TF1* funcB2, Float_t sigmaRangeForSig,Float_t sigmaRangeForBkg, Float_t minMass, Float_t maxMass, Int_t rebin){ |
86d84fa4 | 242 | |
243 | // Compute <pt> from 2D histogram M vs pt | |
244 | ||
245 | //Make 2D histos in the desired pt range | |
246 | Int_t start=hMassD->FindBin(ptmin); | |
247 | Int_t end=hMassD->FindBin(ptmax)-1; | |
248 | const Int_t nx=end-start; | |
249 | TH2F *hMassDpt=new TH2F("hptmass","hptmass",nx,ptmin,ptmax,hMassD->GetNbinsY(),hMassD->GetYaxis()->GetBinLowEdge(1),hMassD->GetYaxis()->GetBinLowEdge(hMassD->GetNbinsY())+hMassD->GetYaxis()->GetBinWidth(hMassD->GetNbinsY())); | |
250 | for(Int_t ix=start;ix<end;ix++){ | |
251 | for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){ | |
252 | hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy)); | |
253 | hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy)); | |
254 | } | |
255 | } | |
256 | ||
257 | Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit; | |
258 | Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit; | |
259 | Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig); | |
260 | Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig); | |
261 | Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig); | |
262 | Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig); | |
263 | // printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin); | |
264 | ||
265 | Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit; | |
19df7ca5 | 266 | Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2); |
86d84fa4 | 267 | Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow); |
268 | Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow); | |
269 | Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow); | |
270 | Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit; | |
271 | Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi); | |
19df7ca5 | 272 | Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1); |
86d84fa4 | 273 | Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi); |
274 | Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi); | |
275 | // printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin); | |
276 | ||
277 | Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin); | |
278 | Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin); | |
279 | Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin); | |
280 | // printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi); | |
281 | ||
282 | TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow); | |
283 | TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi); | |
284 | TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig); | |
285 | ||
286 | hMptBkgLo->Rebin(rebin); | |
287 | hMptBkgHi->Rebin(rebin); | |
288 | hMptSigReg->Rebin(rebin); | |
289 | ||
290 | hMptBkgLo->Sumw2(); | |
291 | hMptBkgHi->Sumw2(); | |
292 | TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin"); | |
293 | hMptBkgLoScal->Scale(bkgSig/bkgLow); | |
294 | TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin"); | |
295 | hMptBkgHiScal->Scale(bkgSig/bkgHi); | |
296 | ||
297 | TH1F* hMptBkgAver=0x0; | |
298 | hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin"); | |
299 | hMptBkgAver->Add(hMptBkgHiScal); | |
300 | hMptBkgAver->Scale(0.5); | |
301 | TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin"); | |
302 | hMptSig->Add(hMptBkgAver,-1.); | |
303 | ||
304 | averagePt = hMptSig->GetMean(); | |
305 | errorPt = hMptSig->GetMeanError(); | |
306 | ||
307 | delete hMptBkgLo; | |
308 | delete hMptBkgHi; | |
309 | delete hMptSigReg; | |
310 | delete hMptBkgLoScal; | |
311 | delete hMptBkgHiScal; | |
312 | delete hMptBkgAver; | |
313 | delete hMassDpt; | |
314 | delete hMptSig; | |
315 | ||
15917f37 | 316 | } |
317 | //____________________________________________________________________________ | |
318 | Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) { | |
319 | // true impact parameter calculation for Dzero | |
320 | ||
321 | if(!partD || !arrayMC || !mcHeader) return 99999.; | |
322 | Int_t code=partD->GetPdgCode(); | |
323 | if(TMath::Abs(code)!=421) return 99999.; | |
324 | ||
325 | Double_t vtxTrue[3]; | |
326 | mcHeader->GetVertex(vtxTrue); | |
327 | Double_t origD[3]; | |
328 | partD->XvYvZv(origD); | |
329 | Short_t charge=partD->Charge(); | |
330 | Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2]; | |
331 | for(Int_t iDau=0; iDau<2; iDau++){ | |
332 | pXdauTrue[iDau]=0.; | |
333 | pYdauTrue[iDau]=0.; | |
334 | pZdauTrue[iDau]=0.; | |
335 | } | |
336 | ||
337 | Int_t nDau=partD->GetNDaughters(); | |
338 | Int_t labelFirstDau = partD->GetDaughter(0); | |
339 | if(nDau==2){ | |
340 | for(Int_t iDau=0; iDau<2; iDau++){ | |
341 | Int_t ind = labelFirstDau+iDau; | |
342 | AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind)); | |
343 | if(!part){ | |
344 | printf("Daughter particle not found in MC array"); | |
345 | return 99999.; | |
346 | } | |
347 | pXdauTrue[iDau]=part->Px(); | |
348 | pYdauTrue[iDau]=part->Py(); | |
349 | pZdauTrue[iDau]=part->Pz(); | |
350 | } | |
351 | }else{ | |
352 | return 99999.; | |
353 | } | |
354 | ||
48efc6ad | 355 | Double_t d0dummy[2]={0.,0.}; |
15917f37 | 356 | AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy); |
357 | return aodDvsMC.ImpParXY(); | |
358 | ||
86d84fa4 | 359 | } |
b176db3b | 360 | //____________________________________________________________________________ |
361 | Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) { | |
362 | // true impact parameter calculation for Dplus | |
363 | ||
364 | if(!partD || !arrayMC || !mcHeader) return 99999.; | |
365 | Int_t code=partD->GetPdgCode(); | |
366 | if(TMath::Abs(code)!=411) return 99999.; | |
367 | ||
368 | Double_t vtxTrue[3]; | |
369 | mcHeader->GetVertex(vtxTrue); | |
370 | Double_t origD[3]; | |
371 | partD->XvYvZv(origD); | |
372 | Short_t charge=partD->Charge(); | |
373 | Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3]; | |
374 | for(Int_t iDau=0; iDau<3; iDau++){ | |
375 | pXdauTrue[iDau]=0.; | |
376 | pYdauTrue[iDau]=0.; | |
377 | pZdauTrue[iDau]=0.; | |
378 | } | |
379 | ||
380 | Int_t nDau=partD->GetNDaughters(); | |
381 | Int_t labelFirstDau = partD->GetDaughter(0); | |
382 | if(nDau==3){ | |
383 | for(Int_t iDau=0; iDau<3; iDau++){ | |
384 | Int_t ind = labelFirstDau+iDau; | |
385 | AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind)); | |
386 | if(!part){ | |
387 | printf("Daughter particle not found in MC array"); | |
388 | return 99999.; | |
389 | } | |
390 | pXdauTrue[iDau]=part->Px(); | |
391 | pYdauTrue[iDau]=part->Py(); | |
392 | pZdauTrue[iDau]=part->Pz(); | |
393 | } | |
394 | }else if(nDau==2){ | |
395 | Int_t theDau=0; | |
396 | for(Int_t iDau=0; iDau<2; iDau++){ | |
397 | Int_t ind = labelFirstDau+iDau; | |
398 | AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind)); | |
399 | if(!part){ | |
400 | printf("Daughter particle not found in MC array"); | |
401 | return 99999.; | |
402 | } | |
403 | Int_t pdgCode=TMath::Abs(part->GetPdgCode()); | |
404 | if(pdgCode==211 || pdgCode==321){ | |
405 | pXdauTrue[theDau]=part->Px(); | |
406 | pYdauTrue[theDau]=part->Py(); | |
407 | pZdauTrue[theDau]=part->Pz(); | |
408 | ++theDau; | |
409 | }else{ | |
410 | Int_t nDauRes=part->GetNDaughters(); | |
411 | if(nDauRes==2){ | |
412 | Int_t labelFirstDauRes = part->GetDaughter(0); | |
413 | for(Int_t iDauRes=0; iDauRes<2; iDauRes++){ | |
414 | Int_t indDR = labelFirstDauRes+iDauRes; | |
415 | AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR)); | |
416 | if(!partDR){ | |
417 | printf("Daughter particle not found in MC array"); | |
418 | return 99999.; | |
419 | } | |
420 | ||
421 | Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode()); | |
422 | if(pdgCodeDR==211 || pdgCodeDR==321){ | |
423 | pXdauTrue[theDau]=partDR->Px(); | |
424 | pYdauTrue[theDau]=partDR->Py(); | |
425 | pZdauTrue[theDau]=partDR->Pz(); | |
426 | ++theDau; | |
427 | } | |
428 | } | |
429 | } | |
430 | } | |
431 | } | |
432 | if(theDau!=3){ | |
433 | printf("Wrong number of decay prongs"); | |
434 | return 99999.; | |
435 | } | |
436 | } | |
437 | ||
438 | Double_t d0dummy[3]={0.,0.,0.}; | |
439 | AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy); | |
440 | return aodDvsMC.ImpParXY(); | |
441 | ||
442 | } | |
443 | ||
444 | ||
445 | ||
686e4710 | 446 | //____________________________________________________________________________ |
447 | Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) { | |
448 | // | |
449 | // Correct the number of accepted tracklets based on a TProfile Hist | |
450 | // | |
451 | // | |
452 | ||
453 | if(TMath::Abs(vtxZ)>10.0){ | |
280cc141 | 454 | // printf("ERROR: Z vertex out of range for correction of multiplicity\n"); |
686e4710 | 455 | return uncorrectedNacc; |
456 | } | |
457 | ||
458 | if(!estimatorAvg){ | |
459 | printf("ERROR: Missing TProfile for correction of multiplicity\n"); | |
460 | return uncorrectedNacc; | |
461 | } | |
462 | ||
463 | Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ)); | |
464 | ||
465 | Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1); | |
466 | ||
467 | Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM)); | |
468 | ||
469 | return correctedNacc; | |
470 | } | |
7f51a627 | 471 | //______________________________________________________________________ |
472 | TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){ | |
f98f4c22 | 473 | // get the name of the generator that produced a given particle |
474 | ||
475 | Int_t nsumpart=0; | |
476 | TList *lh=header->GetCocktailHeaders(); | |
477 | Int_t nh=lh->GetEntries(); | |
478 | for(Int_t i=0;i<nh;i++){ | |
479 | AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i); | |
480 | TString genname=gh->GetName(); | |
481 | Int_t npart=gh->NProduced(); | |
482 | if(label>=nsumpart && label<(nsumpart+npart)) return genname; | |
483 | nsumpart+=npart; | |
484 | } | |
485 | TString empty=""; | |
486 | return empty; | |
487 | } | |
7c0e9154 | 488 | //_____________________________________________________________________ |
489 | void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){ | |
490 | ||
491 | // method to check if a track comes from a given generator | |
f98f4c22 | 492 | |
493 | Int_t lab=TMath::Abs(track->GetLabel()); | |
7c0e9154 | 494 | nameGen=GetGenerator(lab,header); |
f98f4c22 | 495 | |
496 | // Int_t countControl=0; | |
497 | ||
498 | while(nameGen.IsWhitespace()){ | |
499 | AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab); | |
500 | if(!mcpart){ | |
501 | printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab); | |
502 | break; | |
503 | } | |
504 | Int_t mother = mcpart->GetMother(); | |
505 | if(mother<0){ | |
506 | printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n"); | |
507 | break; | |
508 | } | |
509 | lab=mother; | |
510 | nameGen=GetGenerator(mother,header); | |
511 | // countControl++; | |
512 | // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops | |
513 | // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n"); | |
514 | // break; | |
515 | // } | |
516 | } | |
517 | ||
7c0e9154 | 518 | return; |
519 | } | |
520 | //---------------------------------------------------------------------- | |
521 | Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){ | |
522 | // method to check if a track comes from the signal event or from the underlying Hijing event | |
523 | TString nameGen; | |
524 | GetTrackPrimaryGenerator(track,header,arrayMC,nameGen); | |
525 | ||
f98f4c22 | 526 | if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE; |
527 | ||
528 | return kTRUE; | |
529 | } | |
530 | //____________________________________________________________________________ | |
531 | Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){ | |
532 | // method to check if a D meson candidate comes from the signal event or from the underlying Hijing event | |
533 | ||
534 | Int_t nprongs=cand->GetNProngs(); | |
535 | for(Int_t i=0;i<nprongs;i++){ | |
536 | AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i); | |
537 | if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE; | |
538 | } | |
539 | return kFALSE; | |
540 | } | |
c8781624 | 541 | //____________________________________________________________________________ |
d49d50fc | 542 | Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){ |
543 | // method to check if a cascade candidate comes from the signal event or from the underlying Hijing event | |
544 | ||
545 | AliAODTrack* bach = cand->GetBachelor(); | |
546 | if(IsTrackInjected(bach, header, arrayMC)) { | |
547 | AliDebug(2, "Bachelor is injected, the whole candidate is then injected"); | |
548 | return kTRUE; | |
549 | } | |
550 | AliAODv0* v0 = cand->Getv0(); | |
551 | Int_t nprongs = v0->GetNProngs(); | |
552 | for(Int_t i = 0; i < nprongs; i++){ | |
553 | AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i); | |
554 | if(IsTrackInjected(daugh,header,arrayMC)) { | |
555 | AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i)); | |
556 | return kTRUE; | |
557 | } | |
558 | } | |
559 | return kFALSE; | |
560 | } | |
561 | //____________________________________________________________________________ | |
c8781624 | 562 | Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){ |
563 | // checking whether the mother of the particles come from a charm or a bottom quark | |
564 | ||
565 | Int_t pdgGranma = 0; | |
566 | Int_t mother = 0; | |
567 | mother = mcPart->GetMother(); | |
568 | Int_t istep = 0; | |
569 | Int_t abspdgGranma =0; | |
570 | Bool_t isFromB=kFALSE; | |
571 | Bool_t isQuarkFound=kFALSE; | |
572 | while (mother >0 ){ | |
573 | istep++; | |
574 | AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother)); | |
575 | if (mcGranma){ | |
576 | pdgGranma = mcGranma->GetPdgCode(); | |
577 | abspdgGranma = TMath::Abs(pdgGranma); | |
578 | if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){ | |
579 | isFromB=kTRUE; | |
580 | } | |
581 | if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE; | |
582 | mother = mcGranma->GetMother(); | |
583 | }else{ | |
584 | printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!"); | |
585 | break; | |
586 | } | |
587 | } | |
588 | if(searchUpToQuark && !isQuarkFound) return 0; | |
589 | if(isFromB) return 5; | |
590 | else return 4; | |
591 | ||
592 | } | |
593 | //____________________________________________________________________________ | |
594 | Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){ | |
595 | // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases | |
596 | ||
597 | ||
598 | Int_t pdgD=mcPart->GetPdgCode(); | |
599 | if(TMath::Abs(pdgD)!=421) return -1; | |
600 | ||
601 | Int_t nDau=mcPart->GetNDaughters(); | |
602 | ||
603 | if(nDau==2){ | |
604 | Int_t daughter0 = mcPart->GetDaughter(0); | |
605 | Int_t daughter1 = mcPart->GetDaughter(1); | |
606 | AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0)); | |
607 | AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1)); | |
608 | if(!mcPartDaughter0 || !mcPartDaughter1) return -1; | |
609 | arrayDauLab[0]=daughter0; | |
610 | arrayDauLab[1]=daughter1; | |
611 | Int_t pdgCode0=mcPartDaughter0->GetPdgCode(); | |
612 | Int_t pdgCode1=mcPartDaughter1->GetPdgCode(); | |
613 | if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) && | |
614 | !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){ | |
615 | return -1; | |
616 | } | |
617 | if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1; | |
618 | if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1; | |
619 | if((pdgCode0*pdgCode1)>0) return -1; | |
620 | Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px(); | |
621 | Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py(); | |
622 | Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz(); | |
623 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
624 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
625 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
626 | return 1; | |
627 | } | |
628 | ||
629 | if(nDau==3 || nDau==4){ | |
630 | Int_t nKaons=0; | |
631 | Int_t nPions=0; | |
632 | Double_t sumPxDau=0.; | |
633 | Double_t sumPyDau=0.; | |
634 | Double_t sumPzDau=0.; | |
635 | Int_t nFoundKpi=0; | |
636 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
637 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
638 | Int_t indDau = labelFirstDau+iDau; | |
639 | if(indDau<0) return -1; | |
640 | AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau)); | |
641 | if(!dau) return -1; | |
642 | Int_t pdgdau=dau->GetPdgCode(); | |
643 | if(TMath::Abs(pdgdau)==321){ | |
644 | if(pdgD>0 && pdgdau>0) return -1; | |
645 | if(pdgD<0 && pdgdau<0) return -1; | |
646 | nKaons++; | |
647 | sumPxDau+=dau->Px(); | |
648 | sumPyDau+=dau->Py(); | |
649 | sumPzDau+=dau->Pz(); | |
650 | arrayDauLab[nFoundKpi++]=indDau; | |
651 | if(nFoundKpi>4) return -1; | |
652 | }else if(TMath::Abs(pdgdau)==211){ | |
653 | nPions++; | |
654 | sumPxDau+=dau->Px(); | |
655 | sumPyDau+=dau->Py(); | |
656 | sumPzDau+=dau->Pz(); | |
657 | arrayDauLab[nFoundKpi++]=indDau; | |
658 | if(nFoundKpi>4) return -1; | |
659 | }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){ | |
660 | Int_t nResDau=dau->GetNDaughters(); | |
661 | if(nResDau!=2) return -1; | |
662 | Int_t indFirstResDau=dau->GetDaughter(0); | |
663 | for(Int_t resDau=0; resDau<2; resDau++){ | |
664 | Int_t indResDau=indFirstResDau+resDau; | |
665 | if(indResDau<0) return -1; | |
666 | AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau)); | |
667 | if(!resdau) return -1; | |
668 | Int_t pdgresdau=resdau->GetPdgCode(); | |
669 | if(TMath::Abs(pdgresdau)==321){ | |
670 | if(pdgD>0 && pdgresdau>0) return -1; | |
671 | if(pdgD<0 && pdgresdau<0) return -1; | |
672 | nKaons++; | |
673 | sumPxDau+=resdau->Px(); | |
674 | sumPyDau+=resdau->Py(); | |
675 | sumPzDau+=resdau->Pz(); | |
676 | arrayDauLab[nFoundKpi++]=indResDau; | |
677 | if(nFoundKpi>4) return -1; | |
678 | } | |
679 | if(TMath::Abs(pdgresdau)==211){ | |
680 | nPions++; | |
681 | sumPxDau+=resdau->Px(); | |
682 | sumPyDau+=resdau->Py(); | |
683 | sumPzDau+=resdau->Pz(); | |
684 | arrayDauLab[nFoundKpi++]=indResDau; | |
685 | if(nFoundKpi>4) return -1; | |
686 | } | |
687 | } | |
688 | }else{ | |
689 | return -1; | |
690 | } | |
691 | } | |
692 | if(nPions!=3) return -1; | |
693 | if(nKaons!=1) return -1; | |
694 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1; | |
695 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1; | |
696 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1; | |
697 | return 2; | |
698 | } | |
699 | return -1; | |
700 | } | |
701 | //____________________________________________________________________________ | |
702 | Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){ | |
703 | // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases | |
704 | ||
705 | Int_t pdgD=mcPart->GetPdgCode(); | |
706 | if(TMath::Abs(pdgD)!=411) return -1; | |
707 | ||
708 | Int_t nDau=mcPart->GetNDaughters(); | |
709 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
710 | Int_t nKaons=0; | |
711 | Int_t nPions=0; | |
712 | Double_t sumPxDau=0.; | |
713 | Double_t sumPyDau=0.; | |
714 | Double_t sumPzDau=0.; | |
715 | Int_t nFoundKpi=0; | |
716 | ||
717 | if(nDau==3 || nDau==2){ | |
718 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
719 | Int_t indDau = labelFirstDau+iDau; | |
720 | if(indDau<0) return -1; | |
721 | AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau)); | |
722 | if(!dau) return -1; | |
723 | Int_t pdgdau=dau->GetPdgCode(); | |
724 | if(TMath::Abs(pdgdau)==321){ | |
725 | if(pdgD*pdgdau>0) return -1; | |
726 | nKaons++; | |
727 | sumPxDau+=dau->Px(); | |
728 | sumPyDau+=dau->Py(); | |
729 | sumPzDau+=dau->Pz(); | |
730 | arrayDauLab[nFoundKpi++]=indDau; | |
731 | if(nFoundKpi>3) return -1; | |
732 | }else if(TMath::Abs(pdgdau)==211){ | |
733 | if(pdgD*pdgdau<0) return -1; | |
734 | nPions++; | |
735 | sumPxDau+=dau->Px(); | |
736 | sumPyDau+=dau->Py(); | |
737 | sumPzDau+=dau->Pz(); | |
738 | arrayDauLab[nFoundKpi++]=indDau; | |
739 | if(nFoundKpi>3) return -1; | |
740 | }else if(TMath::Abs(pdgdau)==313){ | |
741 | Int_t nResDau=dau->GetNDaughters(); | |
742 | if(nResDau!=2) return -1; | |
743 | Int_t indFirstResDau=dau->GetDaughter(0); | |
744 | for(Int_t resDau=0; resDau<2; resDau++){ | |
745 | Int_t indResDau=indFirstResDau+resDau; | |
746 | if(indResDau<0) return -1; | |
747 | AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau)); | |
748 | if(!resdau) return -1; | |
749 | Int_t pdgresdau=resdau->GetPdgCode(); | |
750 | if(TMath::Abs(pdgresdau)==321){ | |
751 | if(pdgD*pdgresdau>0) return -1; | |
752 | sumPxDau+=resdau->Px(); | |
753 | sumPyDau+=resdau->Py(); | |
754 | sumPzDau+=resdau->Pz(); | |
755 | nKaons++; | |
756 | arrayDauLab[nFoundKpi++]=indResDau; | |
757 | if(nFoundKpi>3) return -1; | |
758 | } | |
759 | if(TMath::Abs(pdgresdau)==211){ | |
760 | if(pdgD*pdgresdau<0) return -1; | |
761 | sumPxDau+=resdau->Px(); | |
762 | sumPyDau+=resdau->Py(); | |
763 | sumPzDau+=resdau->Pz(); | |
764 | nPions++; | |
765 | arrayDauLab[nFoundKpi++]=indResDau; | |
766 | if(nFoundKpi>3) return -1; | |
767 | } | |
768 | } | |
769 | }else{ | |
770 | return -1; | |
771 | } | |
772 | } | |
773 | if(nPions!=2) return -1; | |
774 | if(nKaons!=1) return -1; | |
775 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
776 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
777 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
778 | if(nDau==3) return 1; | |
779 | else if(nDau==2) return 2; | |
780 | } | |
781 | ||
782 | return -1; | |
783 | ||
784 | } |