]>
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> | |
ee9c9d9b | 23 | #include <TParticle.h> |
24 | #include "AliStack.h" | |
a6c5a2e9 | 25 | #include "AliAODEvent.h" |
b176db3b | 26 | #include "AliAODMCHeader.h" |
7f51a627 | 27 | #include "AliGenEventHeader.h" |
25504150 | 28 | #include "AliAODMCParticle.h" |
b176db3b | 29 | #include "AliAODRecoDecayHF.h" |
a6c5a2e9 | 30 | #include "AliVertexingHFUtils.h" |
31 | ||
d947ea9c | 32 | /* $Id$ */ |
a6c5a2e9 | 33 | |
34 | /////////////////////////////////////////////////////////////////// | |
35 | // // | |
36 | // Class with functions useful for different D2H analyses // | |
37 | // - event plane resolution // | |
38 | // - <pt> calculation with side band subtraction // | |
25504150 | 39 | // - tracklet multiplicity calculation // |
a6c5a2e9 | 40 | // Origin: F.Prino, Torino, prino@to.infn.it // |
41 | // // | |
42 | /////////////////////////////////////////////////////////////////// | |
43 | ||
44 | ClassImp(AliVertexingHFUtils) | |
45 | ||
46 | //______________________________________________________________________ | |
47 | AliVertexingHFUtils::AliVertexingHFUtils():TObject(), | |
48 | fK(1), | |
49 | fSubRes(1.), | |
50 | fMinEtaForTracklets(-1.), | |
51 | fMaxEtaForTracklets(1.) | |
52 | { | |
53 | // Default contructor | |
54 | } | |
55 | ||
56 | ||
57 | //______________________________________________________________________ | |
58 | AliVertexingHFUtils::AliVertexingHFUtils(Int_t k): | |
59 | TObject(), | |
60 | fK(k), | |
61 | fSubRes(1.), | |
62 | fMinEtaForTracklets(-1.), | |
63 | fMaxEtaForTracklets(1.) | |
64 | { | |
65 | // Standard constructor | |
66 | } | |
67 | ||
68 | ||
5d1396f9 | 69 | //______________________________________________________________________ |
70 | void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance,Double_t &errsignificance){ | |
71 | // calculate significance from S, B and errors | |
72 | ||
73 | ||
74 | Double_t errSigSq=errsignal*errsignal; | |
75 | Double_t errBkgSq=errbackground*errbackground; | |
76 | Double_t sigPlusBkg=signal+background; | |
77 | if(sigPlusBkg>0. && signal>0.){ | |
78 | significance = signal/TMath::Sqrt(signal+background); | |
79 | errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal); | |
80 | }else{ | |
81 | significance=0.; | |
82 | errsignificance=0.; | |
83 | } | |
84 | return; | |
85 | ||
86 | } | |
a6c5a2e9 | 87 | //______________________________________________________________________ |
88 | Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){ | |
89 | // compute chi from polynomial approximation | |
90 | Double_t c[5]; | |
91 | if(k==1){ | |
92 | c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283; | |
93 | } | |
94 | else if(k==2){ | |
95 | c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815; | |
96 | } | |
97 | else return -1; | |
98 | 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; | |
99 | } | |
100 | ||
101 | //______________________________________________________________________ | |
102 | Double_t AliVertexingHFUtils:: ResolK1(Double_t x){ | |
103 | return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4)); | |
104 | } | |
105 | ||
106 | ||
107 | //______________________________________________________________________ | |
108 | Double_t AliVertexingHFUtils::FindChi(Double_t res, Int_t k){ | |
109 | // compute chi variable (=v2*sqrt(N)) from external values | |
110 | ||
111 | Double_t x1=0; | |
112 | Double_t x2=15; | |
113 | Double_t y1,y2; | |
114 | if(k==1){ | |
115 | y1=ResolK1(x1)-res; | |
116 | y2=ResolK1(x2)-res; | |
117 | } | |
118 | else if(k==2){ | |
119 | y1=Pol(x1,2)-res; | |
120 | y2=Pol(x2,2)-res; | |
121 | } | |
122 | else return -1; | |
123 | ||
124 | if(y1*y2>0) return -1; | |
125 | if(y1==0) return y1; | |
126 | if(y2==0) return y2; | |
127 | Double_t xmed,ymed; | |
128 | Int_t jiter=0; | |
129 | while((x2-x1)>0.0001){ | |
130 | xmed=0.5*(x1+x2); | |
131 | if(k==1){ | |
132 | y1=ResolK1(x1)-res; | |
133 | y2=ResolK1(x2)-res; | |
134 | ymed=ResolK1(xmed)-res; | |
135 | } | |
136 | else if(k==2){ | |
137 | y1=Pol(x1,2)-res; | |
138 | y2=Pol(x2,2)-res; | |
139 | ymed=Pol(xmed,2)-res; | |
140 | } | |
141 | else return -1; | |
142 | if((y1*ymed)<0) x2=xmed; | |
143 | if((y2*ymed)<0)x1=xmed; | |
144 | if(ymed==0) return xmed; | |
145 | jiter++; | |
146 | } | |
147 | return 0.5*(x1+x2); | |
148 | } | |
149 | ||
150 | //______________________________________________________________________ | |
151 | Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){ | |
152 | // computes event plane resolution starting from sub event resolution | |
153 | Double_t chisub=FindChi(resSub,k); | |
154 | Double_t chifull=chisub*TMath::Sqrt(2); | |
155 | if(k==1) return ResolK1(chifull); | |
156 | else if(k==2) return Pol(chifull,2); | |
157 | else{ | |
158 | printf("k should be <=2\n"); | |
159 | return 1.; | |
160 | } | |
161 | } | |
162 | ||
163 | //______________________________________________________________________ | |
164 | Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){ | |
165 | // computes event plane resolution starting from sub event correlation histogram | |
166 | if(!hSubEvCorr) return 1.; | |
167 | Double_t resSub=GetSubEvResol(hSubEvCorr); | |
168 | return GetFullEvResol(resSub,k); | |
169 | } | |
170 | //______________________________________________________________________ | |
171 | Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){ | |
172 | // computes low limit event plane resolution starting from sub event correlation histogram | |
173 | if(!hSubEvCorr) return 1.; | |
174 | Double_t resSub=GetSubEvResolLowLim(hSubEvCorr); | |
175 | printf("%f\n",resSub); | |
176 | return GetFullEvResol(resSub,k); | |
177 | } | |
178 | //______________________________________________________________________ | |
179 | Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){ | |
180 | // computes high limit event plane resolution starting from sub event correlation histogram | |
181 | if(!hSubEvCorr) return 1.; | |
182 | Double_t resSub=GetSubEvResolHighLim(hSubEvCorr); | |
183 | printf("%f\n",resSub); | |
184 | return GetFullEvResol(resSub,k); | |
185 | } | |
186 | //______________________________________________________________________ | |
187 | Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){ | |
188 | // counts tracklets in given eta range | |
189 | AliAODTracklets* tracklets=ev->GetTracklets(); | |
190 | Int_t nTr=tracklets->GetNumberOfTracklets(); | |
191 | Int_t count=0; | |
192 | for(Int_t iTr=0; iTr<nTr; iTr++){ | |
193 | Double_t theta=tracklets->GetTheta(iTr); | |
194 | Double_t eta=-TMath::Log(TMath::Tan(theta/2.)); | |
195 | if(eta>mineta && eta<maxeta) count++; | |
196 | } | |
197 | return count; | |
198 | } | |
86d84fa4 | 199 | //______________________________________________________________________ |
25504150 | 200 | Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){ |
201 | // counts generated particles in fgiven eta range | |
202 | ||
203 | Int_t nChargedMC=0; | |
204 | for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){ | |
205 | AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i); | |
206 | Int_t charge = part->Charge(); | |
207 | Double_t eta = part->Eta(); | |
208 | if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++; | |
209 | } | |
210 | return nChargedMC; | |
211 | } | |
212 | //______________________________________________________________________ | |
213 | Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){ | |
214 | // counts generated primary particles in given eta range | |
215 | ||
216 | Int_t nChargedMC=0; | |
217 | for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){ | |
218 | AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i); | |
219 | Int_t charge = part->Charge(); | |
220 | Double_t eta = part->Eta(); | |
221 | if(charge!=0 && eta>mineta && eta<maxeta){ | |
222 | if(part->IsPrimary())nChargedMC++; | |
223 | } | |
224 | } | |
225 | return nChargedMC; | |
226 | } | |
227 | //______________________________________________________________________ | |
228 | Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){ | |
229 | // counts generated primary particles in given eta range | |
230 | ||
231 | Int_t nChargedMC=0; | |
232 | for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){ | |
233 | AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i); | |
234 | Int_t charge = part->Charge(); | |
235 | Double_t eta = part->Eta(); | |
236 | if(charge!=0 && eta>mineta && eta<maxeta){ | |
237 | if(part->IsPhysicalPrimary())nChargedMC++; | |
238 | } | |
239 | } | |
240 | return nChargedMC; | |
241 | } | |
b242bddf | 242 | |
243 | ||
244 | //______________________________________________________________________ | |
245 | Double_t AliVertexingHFUtils::GetVZEROAEqualizedMultiplicity(AliAODEvent* ev){ | |
246 | // | |
247 | // Method to get VZERO-A equalized multiplicty as done in AliCentralitySelectionTask | |
248 | // getting the equalized VZERO factors from the tender or AOD | |
249 | // http://git.cern.ch/pubweb/AliRoot.git/blob/HEAD:/ANALYSIS/AliCentralitySelectionTask.cxx#l1345 | |
250 | ||
251 | Double_t multV0AEq=0; | |
252 | for(Int_t iCh = 32; iCh < 64; ++iCh) { | |
253 | Double_t mult = ev->GetVZEROEqMultiplicity(iCh); | |
254 | multV0AEq += mult; | |
255 | } | |
256 | return multV0AEq; | |
257 | } | |
258 | ||
259 | //______________________________________________________________________ | |
260 | Double_t AliVertexingHFUtils::GetVZEROCEqualizedMultiplicity(AliAODEvent* ev){ | |
261 | // Method to get VZERO-C equalized multiplicty as done in AliCentralitySelectionTask | |
262 | // getting the equalized VZERO factors from the tender or AOD | |
263 | // http://git.cern.ch/pubweb/AliRoot.git/blob/HEAD:/ANALYSIS/AliCentralitySelectionTask.cxx#l1345 | |
264 | ||
265 | Double_t multV0CEq=0; | |
266 | for(Int_t iCh = 0; iCh < 32; ++iCh) { | |
267 | Double_t mult = ev->GetVZEROEqMultiplicity(iCh); | |
268 | multV0CEq += mult; | |
269 | } | |
270 | return multV0CEq; | |
271 | } | |
272 | ||
25504150 | 273 | //______________________________________________________________________ |
19df7ca5 | 274 | 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 | 275 | |
276 | // Compute <pt> from 2D histogram M vs pt | |
277 | ||
278 | //Make 2D histos in the desired pt range | |
279 | Int_t start=hMassD->FindBin(ptmin); | |
280 | Int_t end=hMassD->FindBin(ptmax)-1; | |
281 | const Int_t nx=end-start; | |
282 | 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())); | |
283 | for(Int_t ix=start;ix<end;ix++){ | |
284 | for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){ | |
285 | hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy)); | |
286 | hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy)); | |
287 | } | |
288 | } | |
289 | ||
290 | Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit; | |
291 | Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit; | |
292 | Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig); | |
293 | Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig); | |
294 | Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig); | |
295 | Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig); | |
296 | // printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin); | |
297 | ||
298 | Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit; | |
19df7ca5 | 299 | Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2); |
86d84fa4 | 300 | Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow); |
301 | Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow); | |
302 | Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow); | |
303 | Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit; | |
304 | Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi); | |
19df7ca5 | 305 | Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1); |
86d84fa4 | 306 | Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi); |
307 | Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi); | |
308 | // printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin); | |
309 | ||
310 | Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin); | |
311 | Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin); | |
312 | Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin); | |
313 | // printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi); | |
314 | ||
315 | TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow); | |
316 | TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi); | |
317 | TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig); | |
318 | ||
319 | hMptBkgLo->Rebin(rebin); | |
320 | hMptBkgHi->Rebin(rebin); | |
321 | hMptSigReg->Rebin(rebin); | |
322 | ||
323 | hMptBkgLo->Sumw2(); | |
324 | hMptBkgHi->Sumw2(); | |
325 | TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin"); | |
326 | hMptBkgLoScal->Scale(bkgSig/bkgLow); | |
327 | TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin"); | |
328 | hMptBkgHiScal->Scale(bkgSig/bkgHi); | |
329 | ||
330 | TH1F* hMptBkgAver=0x0; | |
331 | hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin"); | |
332 | hMptBkgAver->Add(hMptBkgHiScal); | |
333 | hMptBkgAver->Scale(0.5); | |
334 | TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin"); | |
335 | hMptSig->Add(hMptBkgAver,-1.); | |
336 | ||
337 | averagePt = hMptSig->GetMean(); | |
338 | errorPt = hMptSig->GetMeanError(); | |
339 | ||
340 | delete hMptBkgLo; | |
341 | delete hMptBkgHi; | |
342 | delete hMptSigReg; | |
343 | delete hMptBkgLoScal; | |
344 | delete hMptBkgHiScal; | |
345 | delete hMptBkgAver; | |
346 | delete hMassDpt; | |
347 | delete hMptSig; | |
348 | ||
15917f37 | 349 | } |
350 | //____________________________________________________________________________ | |
e6eaae4c | 351 | Bool_t AliVertexingHFUtils::CheckT0TriggerFired(AliAODEvent* aodEv){ |
352 | // check if T0VTX trigger was fired, based on a workaround suggested by Alla | |
353 | const Double32_t *mean = aodEv->GetT0TOF(); | |
94d7319f | 354 | if(mean && mean[0]<9999.) return kTRUE; |
e6eaae4c | 355 | else return kFALSE; |
356 | } | |
357 | //____________________________________________________________________________ | |
15917f37 | 358 | Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) { |
359 | // true impact parameter calculation for Dzero | |
360 | ||
361 | if(!partD || !arrayMC || !mcHeader) return 99999.; | |
362 | Int_t code=partD->GetPdgCode(); | |
363 | if(TMath::Abs(code)!=421) return 99999.; | |
364 | ||
365 | Double_t vtxTrue[3]; | |
366 | mcHeader->GetVertex(vtxTrue); | |
367 | Double_t origD[3]; | |
368 | partD->XvYvZv(origD); | |
369 | Short_t charge=partD->Charge(); | |
370 | Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2]; | |
371 | for(Int_t iDau=0; iDau<2; iDau++){ | |
372 | pXdauTrue[iDau]=0.; | |
373 | pYdauTrue[iDau]=0.; | |
374 | pZdauTrue[iDau]=0.; | |
375 | } | |
376 | ||
377 | Int_t nDau=partD->GetNDaughters(); | |
378 | Int_t labelFirstDau = partD->GetDaughter(0); | |
379 | if(nDau==2){ | |
380 | for(Int_t iDau=0; iDau<2; iDau++){ | |
381 | Int_t ind = labelFirstDau+iDau; | |
382 | AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind)); | |
383 | if(!part){ | |
384 | printf("Daughter particle not found in MC array"); | |
385 | return 99999.; | |
386 | } | |
387 | pXdauTrue[iDau]=part->Px(); | |
388 | pYdauTrue[iDau]=part->Py(); | |
389 | pZdauTrue[iDau]=part->Pz(); | |
390 | } | |
391 | }else{ | |
392 | return 99999.; | |
393 | } | |
394 | ||
48efc6ad | 395 | Double_t d0dummy[2]={0.,0.}; |
15917f37 | 396 | AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy); |
397 | return aodDvsMC.ImpParXY(); | |
398 | ||
86d84fa4 | 399 | } |
b176db3b | 400 | //____________________________________________________________________________ |
401 | Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) { | |
402 | // true impact parameter calculation for Dplus | |
403 | ||
404 | if(!partD || !arrayMC || !mcHeader) return 99999.; | |
405 | Int_t code=partD->GetPdgCode(); | |
406 | if(TMath::Abs(code)!=411) return 99999.; | |
407 | ||
408 | Double_t vtxTrue[3]; | |
409 | mcHeader->GetVertex(vtxTrue); | |
410 | Double_t origD[3]; | |
411 | partD->XvYvZv(origD); | |
412 | Short_t charge=partD->Charge(); | |
413 | Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3]; | |
414 | for(Int_t iDau=0; iDau<3; iDau++){ | |
415 | pXdauTrue[iDau]=0.; | |
416 | pYdauTrue[iDau]=0.; | |
417 | pZdauTrue[iDau]=0.; | |
418 | } | |
419 | ||
420 | Int_t nDau=partD->GetNDaughters(); | |
421 | Int_t labelFirstDau = partD->GetDaughter(0); | |
422 | if(nDau==3){ | |
423 | for(Int_t iDau=0; iDau<3; iDau++){ | |
424 | Int_t ind = labelFirstDau+iDau; | |
425 | AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind)); | |
426 | if(!part){ | |
427 | printf("Daughter particle not found in MC array"); | |
428 | return 99999.; | |
429 | } | |
430 | pXdauTrue[iDau]=part->Px(); | |
431 | pYdauTrue[iDau]=part->Py(); | |
432 | pZdauTrue[iDau]=part->Pz(); | |
433 | } | |
434 | }else if(nDau==2){ | |
435 | Int_t theDau=0; | |
436 | for(Int_t iDau=0; iDau<2; iDau++){ | |
437 | Int_t ind = labelFirstDau+iDau; | |
438 | AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind)); | |
439 | if(!part){ | |
440 | printf("Daughter particle not found in MC array"); | |
441 | return 99999.; | |
442 | } | |
443 | Int_t pdgCode=TMath::Abs(part->GetPdgCode()); | |
444 | if(pdgCode==211 || pdgCode==321){ | |
445 | pXdauTrue[theDau]=part->Px(); | |
446 | pYdauTrue[theDau]=part->Py(); | |
447 | pZdauTrue[theDau]=part->Pz(); | |
448 | ++theDau; | |
449 | }else{ | |
450 | Int_t nDauRes=part->GetNDaughters(); | |
451 | if(nDauRes==2){ | |
452 | Int_t labelFirstDauRes = part->GetDaughter(0); | |
453 | for(Int_t iDauRes=0; iDauRes<2; iDauRes++){ | |
454 | Int_t indDR = labelFirstDauRes+iDauRes; | |
455 | AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR)); | |
456 | if(!partDR){ | |
457 | printf("Daughter particle not found in MC array"); | |
458 | return 99999.; | |
459 | } | |
460 | ||
461 | Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode()); | |
462 | if(pdgCodeDR==211 || pdgCodeDR==321){ | |
463 | pXdauTrue[theDau]=partDR->Px(); | |
464 | pYdauTrue[theDau]=partDR->Py(); | |
465 | pZdauTrue[theDau]=partDR->Pz(); | |
466 | ++theDau; | |
467 | } | |
468 | } | |
469 | } | |
470 | } | |
471 | } | |
472 | if(theDau!=3){ | |
473 | printf("Wrong number of decay prongs"); | |
474 | return 99999.; | |
475 | } | |
476 | } | |
477 | ||
478 | Double_t d0dummy[3]={0.,0.,0.}; | |
479 | AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy); | |
480 | return aodDvsMC.ImpParXY(); | |
481 | ||
482 | } | |
483 | ||
484 | ||
485 | ||
686e4710 | 486 | //____________________________________________________________________________ |
487 | Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) { | |
488 | // | |
489 | // Correct the number of accepted tracklets based on a TProfile Hist | |
490 | // | |
491 | // | |
492 | ||
493 | if(TMath::Abs(vtxZ)>10.0){ | |
280cc141 | 494 | // printf("ERROR: Z vertex out of range for correction of multiplicity\n"); |
686e4710 | 495 | return uncorrectedNacc; |
496 | } | |
497 | ||
498 | if(!estimatorAvg){ | |
499 | printf("ERROR: Missing TProfile for correction of multiplicity\n"); | |
500 | return uncorrectedNacc; | |
501 | } | |
502 | ||
503 | Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ)); | |
504 | ||
505 | Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1); | |
506 | ||
507 | Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM)); | |
508 | ||
509 | return correctedNacc; | |
510 | } | |
7f51a627 | 511 | //______________________________________________________________________ |
512 | TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){ | |
f98f4c22 | 513 | // get the name of the generator that produced a given particle |
514 | ||
515 | Int_t nsumpart=0; | |
516 | TList *lh=header->GetCocktailHeaders(); | |
517 | Int_t nh=lh->GetEntries(); | |
518 | for(Int_t i=0;i<nh;i++){ | |
519 | AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i); | |
520 | TString genname=gh->GetName(); | |
521 | Int_t npart=gh->NProduced(); | |
522 | if(label>=nsumpart && label<(nsumpart+npart)) return genname; | |
523 | nsumpart+=npart; | |
524 | } | |
525 | TString empty=""; | |
526 | return empty; | |
527 | } | |
7c0e9154 | 528 | //_____________________________________________________________________ |
529 | void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){ | |
530 | ||
531 | // method to check if a track comes from a given generator | |
f98f4c22 | 532 | |
533 | Int_t lab=TMath::Abs(track->GetLabel()); | |
7c0e9154 | 534 | nameGen=GetGenerator(lab,header); |
f98f4c22 | 535 | |
536 | // Int_t countControl=0; | |
537 | ||
538 | while(nameGen.IsWhitespace()){ | |
539 | AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab); | |
540 | if(!mcpart){ | |
541 | printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab); | |
542 | break; | |
543 | } | |
544 | Int_t mother = mcpart->GetMother(); | |
545 | if(mother<0){ | |
546 | printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n"); | |
547 | break; | |
548 | } | |
549 | lab=mother; | |
550 | nameGen=GetGenerator(mother,header); | |
551 | // countControl++; | |
552 | // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops | |
553 | // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n"); | |
554 | // break; | |
555 | // } | |
556 | } | |
557 | ||
7c0e9154 | 558 | return; |
559 | } | |
560 | //---------------------------------------------------------------------- | |
561 | Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){ | |
562 | // method to check if a track comes from the signal event or from the underlying Hijing event | |
563 | TString nameGen; | |
564 | GetTrackPrimaryGenerator(track,header,arrayMC,nameGen); | |
565 | ||
f98f4c22 | 566 | if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE; |
567 | ||
568 | return kTRUE; | |
569 | } | |
570 | //____________________________________________________________________________ | |
571 | Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){ | |
572 | // method to check if a D meson candidate comes from the signal event or from the underlying Hijing event | |
573 | ||
574 | Int_t nprongs=cand->GetNProngs(); | |
575 | for(Int_t i=0;i<nprongs;i++){ | |
576 | AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i); | |
577 | if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE; | |
578 | } | |
579 | return kFALSE; | |
580 | } | |
c8781624 | 581 | //____________________________________________________________________________ |
d49d50fc | 582 | Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){ |
583 | // method to check if a cascade candidate comes from the signal event or from the underlying Hijing event | |
584 | ||
585 | AliAODTrack* bach = cand->GetBachelor(); | |
586 | if(IsTrackInjected(bach, header, arrayMC)) { | |
587 | AliDebug(2, "Bachelor is injected, the whole candidate is then injected"); | |
588 | return kTRUE; | |
589 | } | |
590 | AliAODv0* v0 = cand->Getv0(); | |
591 | Int_t nprongs = v0->GetNProngs(); | |
592 | for(Int_t i = 0; i < nprongs; i++){ | |
593 | AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i); | |
594 | if(IsTrackInjected(daugh,header,arrayMC)) { | |
595 | AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i)); | |
596 | return kTRUE; | |
597 | } | |
598 | } | |
599 | return kFALSE; | |
600 | } | |
601 | //____________________________________________________________________________ | |
ee9c9d9b | 602 | Int_t AliVertexingHFUtils::CheckOrigin(AliStack* stack, TParticle *mcPart, Bool_t searchUpToQuark){ |
603 | // checking whether the mother of the particles come from a charm or a bottom quark | |
604 | ||
605 | Int_t pdgGranma = 0; | |
606 | Int_t mother = 0; | |
607 | mother = mcPart->GetFirstMother(); | |
608 | Int_t istep = 0; | |
609 | Int_t abspdgGranma =0; | |
610 | Bool_t isFromB=kFALSE; | |
611 | Bool_t isQuarkFound=kFALSE; | |
612 | while (mother >0 ){ | |
613 | istep++; | |
614 | TParticle* mcGranma = stack->Particle(mother); | |
615 | if (mcGranma){ | |
616 | pdgGranma = mcGranma->GetPdgCode(); | |
617 | abspdgGranma = TMath::Abs(pdgGranma); | |
618 | if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){ | |
619 | isFromB=kTRUE; | |
620 | } | |
621 | if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE; | |
622 | mother = mcGranma->GetFirstMother(); | |
623 | }else{ | |
624 | printf("CheckOrigin: Failed casting the mother particle!"); | |
625 | break; | |
626 | } | |
627 | } | |
628 | if(searchUpToQuark && !isQuarkFound) return 0; | |
629 | if(isFromB) return 5; | |
630 | else return 4; | |
631 | ||
632 | } | |
633 | //____________________________________________________________________________ | |
c8781624 | 634 | Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){ |
635 | // checking whether the mother of the particles come from a charm or a bottom quark | |
636 | ||
637 | Int_t pdgGranma = 0; | |
638 | Int_t mother = 0; | |
639 | mother = mcPart->GetMother(); | |
640 | Int_t istep = 0; | |
641 | Int_t abspdgGranma =0; | |
642 | Bool_t isFromB=kFALSE; | |
643 | Bool_t isQuarkFound=kFALSE; | |
644 | while (mother >0 ){ | |
645 | istep++; | |
646 | AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother)); | |
647 | if (mcGranma){ | |
648 | pdgGranma = mcGranma->GetPdgCode(); | |
649 | abspdgGranma = TMath::Abs(pdgGranma); | |
650 | if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){ | |
651 | isFromB=kTRUE; | |
652 | } | |
653 | if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE; | |
654 | mother = mcGranma->GetMother(); | |
655 | }else{ | |
656 | printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!"); | |
657 | break; | |
658 | } | |
659 | } | |
660 | if(searchUpToQuark && !isQuarkFound) return 0; | |
661 | if(isFromB) return 5; | |
662 | else return 4; | |
663 | ||
664 | } | |
ee9c9d9b | 665 | |
c8781624 | 666 | //____________________________________________________________________________ |
ee9c9d9b | 667 | Int_t AliVertexingHFUtils::CheckD0Decay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ |
c8781624 | 668 | // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases |
669 | ||
ee9c9d9b | 670 | if(label<0) return -1; |
671 | TParticle* mcPart = stack->Particle(label); | |
672 | Int_t pdgD=mcPart->GetPdgCode(); | |
673 | if(TMath::Abs(pdgD)!=421) return -1; | |
674 | ||
675 | Int_t nDau=mcPart->GetNDaughters(); | |
676 | ||
677 | if(nDau==2){ | |
678 | Int_t daughter0 = mcPart->GetDaughter(0); | |
679 | Int_t daughter1 = mcPart->GetDaughter(1); | |
680 | TParticle* mcPartDaughter0 = stack->Particle(daughter0); | |
681 | TParticle* mcPartDaughter1 = stack->Particle(daughter1); | |
682 | if(!mcPartDaughter0 || !mcPartDaughter1) return -1; | |
683 | arrayDauLab[0]=daughter0; | |
684 | arrayDauLab[1]=daughter1; | |
685 | Int_t pdgCode0=mcPartDaughter0->GetPdgCode(); | |
686 | Int_t pdgCode1=mcPartDaughter1->GetPdgCode(); | |
687 | if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) && | |
688 | !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){ | |
689 | return -1; | |
690 | } | |
691 | if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1; | |
692 | if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1; | |
693 | if((pdgCode0*pdgCode1)>0) return -1; | |
694 | Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px(); | |
695 | Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py(); | |
696 | Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz(); | |
697 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
698 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
699 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
700 | return 1; | |
701 | } | |
702 | ||
703 | if(nDau==3 || nDau==4){ | |
704 | Int_t nKaons=0; | |
705 | Int_t nPions=0; | |
706 | Double_t sumPxDau=0.; | |
707 | Double_t sumPyDau=0.; | |
708 | Double_t sumPzDau=0.; | |
709 | Int_t nFoundKpi=0; | |
710 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
711 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
712 | Int_t indDau = labelFirstDau+iDau; | |
713 | if(indDau<0) return -1; | |
714 | TParticle* dau=stack->Particle(indDau); | |
715 | if(!dau) return -1; | |
716 | Int_t pdgdau=dau->GetPdgCode(); | |
717 | if(TMath::Abs(pdgdau)==321){ | |
718 | if(pdgD>0 && pdgdau>0) return -1; | |
719 | if(pdgD<0 && pdgdau<0) return -1; | |
720 | nKaons++; | |
721 | sumPxDau+=dau->Px(); | |
722 | sumPyDau+=dau->Py(); | |
723 | sumPzDau+=dau->Pz(); | |
724 | arrayDauLab[nFoundKpi++]=indDau; | |
725 | if(nFoundKpi>4) return -1; | |
726 | }else if(TMath::Abs(pdgdau)==211){ | |
727 | nPions++; | |
728 | sumPxDau+=dau->Px(); | |
729 | sumPyDau+=dau->Py(); | |
730 | sumPzDau+=dau->Pz(); | |
731 | arrayDauLab[nFoundKpi++]=indDau; | |
732 | if(nFoundKpi>4) return -1; | |
733 | }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){ | |
734 | Int_t nResDau=dau->GetNDaughters(); | |
735 | if(nResDau!=2) return -1; | |
736 | Int_t indFirstResDau=dau->GetDaughter(0); | |
737 | for(Int_t resDau=0; resDau<2; resDau++){ | |
738 | Int_t indResDau=indFirstResDau+resDau; | |
739 | if(indResDau<0) return -1; | |
740 | TParticle* resdau=stack->Particle(indResDau); | |
741 | if(!resdau) return -1; | |
742 | Int_t pdgresdau=resdau->GetPdgCode(); | |
743 | if(TMath::Abs(pdgresdau)==321){ | |
744 | if(pdgD>0 && pdgresdau>0) return -1; | |
745 | if(pdgD<0 && pdgresdau<0) return -1; | |
746 | nKaons++; | |
747 | sumPxDau+=resdau->Px(); | |
748 | sumPyDau+=resdau->Py(); | |
749 | sumPzDau+=resdau->Pz(); | |
750 | arrayDauLab[nFoundKpi++]=indResDau; | |
751 | if(nFoundKpi>4) return -1; | |
752 | } | |
753 | if(TMath::Abs(pdgresdau)==211){ | |
754 | nPions++; | |
755 | sumPxDau+=resdau->Px(); | |
756 | sumPyDau+=resdau->Py(); | |
757 | sumPzDau+=resdau->Pz(); | |
758 | arrayDauLab[nFoundKpi++]=indResDau; | |
759 | if(nFoundKpi>4) return -1; | |
760 | } | |
761 | } | |
762 | }else{ | |
763 | return -1; | |
764 | } | |
765 | } | |
766 | if(nPions!=3) return -1; | |
767 | if(nKaons!=1) return -1; | |
768 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1; | |
769 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1; | |
770 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1; | |
771 | return 2; | |
772 | } | |
773 | return -1; | |
774 | } | |
775 | ||
776 | //____________________________________________________________________________ | |
777 | Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){ | |
778 | // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases | |
c8781624 | 779 | |
780 | Int_t pdgD=mcPart->GetPdgCode(); | |
781 | if(TMath::Abs(pdgD)!=421) return -1; | |
782 | ||
783 | Int_t nDau=mcPart->GetNDaughters(); | |
784 | ||
785 | if(nDau==2){ | |
786 | Int_t daughter0 = mcPart->GetDaughter(0); | |
787 | Int_t daughter1 = mcPart->GetDaughter(1); | |
788 | AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0)); | |
789 | AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1)); | |
790 | if(!mcPartDaughter0 || !mcPartDaughter1) return -1; | |
791 | arrayDauLab[0]=daughter0; | |
792 | arrayDauLab[1]=daughter1; | |
793 | Int_t pdgCode0=mcPartDaughter0->GetPdgCode(); | |
794 | Int_t pdgCode1=mcPartDaughter1->GetPdgCode(); | |
795 | if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) && | |
796 | !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){ | |
797 | return -1; | |
798 | } | |
799 | if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1; | |
800 | if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1; | |
801 | if((pdgCode0*pdgCode1)>0) return -1; | |
802 | Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px(); | |
803 | Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py(); | |
804 | Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz(); | |
805 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
806 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
807 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
808 | return 1; | |
809 | } | |
810 | ||
811 | if(nDau==3 || nDau==4){ | |
812 | Int_t nKaons=0; | |
813 | Int_t nPions=0; | |
814 | Double_t sumPxDau=0.; | |
815 | Double_t sumPyDau=0.; | |
816 | Double_t sumPzDau=0.; | |
817 | Int_t nFoundKpi=0; | |
818 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
819 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
820 | Int_t indDau = labelFirstDau+iDau; | |
821 | if(indDau<0) return -1; | |
822 | AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau)); | |
823 | if(!dau) return -1; | |
824 | Int_t pdgdau=dau->GetPdgCode(); | |
825 | if(TMath::Abs(pdgdau)==321){ | |
826 | if(pdgD>0 && pdgdau>0) return -1; | |
827 | if(pdgD<0 && pdgdau<0) return -1; | |
828 | nKaons++; | |
829 | sumPxDau+=dau->Px(); | |
830 | sumPyDau+=dau->Py(); | |
831 | sumPzDau+=dau->Pz(); | |
832 | arrayDauLab[nFoundKpi++]=indDau; | |
833 | if(nFoundKpi>4) return -1; | |
834 | }else if(TMath::Abs(pdgdau)==211){ | |
835 | nPions++; | |
836 | sumPxDau+=dau->Px(); | |
837 | sumPyDau+=dau->Py(); | |
838 | sumPzDau+=dau->Pz(); | |
839 | arrayDauLab[nFoundKpi++]=indDau; | |
840 | if(nFoundKpi>4) return -1; | |
841 | }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){ | |
842 | Int_t nResDau=dau->GetNDaughters(); | |
843 | if(nResDau!=2) return -1; | |
844 | Int_t indFirstResDau=dau->GetDaughter(0); | |
845 | for(Int_t resDau=0; resDau<2; resDau++){ | |
846 | Int_t indResDau=indFirstResDau+resDau; | |
847 | if(indResDau<0) return -1; | |
848 | AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau)); | |
849 | if(!resdau) return -1; | |
850 | Int_t pdgresdau=resdau->GetPdgCode(); | |
851 | if(TMath::Abs(pdgresdau)==321){ | |
852 | if(pdgD>0 && pdgresdau>0) return -1; | |
853 | if(pdgD<0 && pdgresdau<0) return -1; | |
854 | nKaons++; | |
855 | sumPxDau+=resdau->Px(); | |
856 | sumPyDau+=resdau->Py(); | |
857 | sumPzDau+=resdau->Pz(); | |
858 | arrayDauLab[nFoundKpi++]=indResDau; | |
859 | if(nFoundKpi>4) return -1; | |
860 | } | |
861 | if(TMath::Abs(pdgresdau)==211){ | |
862 | nPions++; | |
863 | sumPxDau+=resdau->Px(); | |
864 | sumPyDau+=resdau->Py(); | |
865 | sumPzDau+=resdau->Pz(); | |
866 | arrayDauLab[nFoundKpi++]=indResDau; | |
867 | if(nFoundKpi>4) return -1; | |
868 | } | |
869 | } | |
870 | }else{ | |
871 | return -1; | |
872 | } | |
873 | } | |
874 | if(nPions!=3) return -1; | |
875 | if(nKaons!=1) return -1; | |
876 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1; | |
877 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1; | |
878 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1; | |
879 | return 2; | |
880 | } | |
881 | return -1; | |
882 | } | |
ee9c9d9b | 883 | //____________________________________________________________________________ |
884 | Int_t AliVertexingHFUtils::CheckDplusDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ | |
885 | // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases | |
886 | ||
887 | if(label<0) return -1; | |
888 | TParticle* mcPart = stack->Particle(label); | |
889 | Int_t pdgD=mcPart->GetPdgCode(); | |
890 | if(TMath::Abs(pdgD)!=411) return -1; | |
891 | ||
892 | Int_t nDau=mcPart->GetNDaughters(); | |
893 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
894 | Int_t nKaons=0; | |
895 | Int_t nPions=0; | |
896 | Double_t sumPxDau=0.; | |
897 | Double_t sumPyDau=0.; | |
898 | Double_t sumPzDau=0.; | |
899 | Int_t nFoundKpi=0; | |
900 | ||
901 | if(nDau==3 || nDau==2){ | |
902 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
903 | Int_t indDau = labelFirstDau+iDau; | |
904 | if(indDau<0) return -1; | |
905 | TParticle* dau=stack->Particle(indDau); | |
906 | if(!dau) return -1; | |
907 | Int_t pdgdau=dau->GetPdgCode(); | |
908 | if(TMath::Abs(pdgdau)==321){ | |
909 | if(pdgD*pdgdau>0) return -1; | |
910 | nKaons++; | |
911 | sumPxDau+=dau->Px(); | |
912 | sumPyDau+=dau->Py(); | |
913 | sumPzDau+=dau->Pz(); | |
914 | arrayDauLab[nFoundKpi++]=indDau; | |
915 | if(nFoundKpi>3) return -1; | |
916 | }else if(TMath::Abs(pdgdau)==211){ | |
917 | if(pdgD*pdgdau<0) return -1; | |
918 | nPions++; | |
919 | sumPxDau+=dau->Px(); | |
920 | sumPyDau+=dau->Py(); | |
921 | sumPzDau+=dau->Pz(); | |
922 | arrayDauLab[nFoundKpi++]=indDau; | |
923 | if(nFoundKpi>3) return -1; | |
924 | }else if(TMath::Abs(pdgdau)==313){ | |
925 | Int_t nResDau=dau->GetNDaughters(); | |
926 | if(nResDau!=2) return -1; | |
927 | Int_t indFirstResDau=dau->GetDaughter(0); | |
928 | for(Int_t resDau=0; resDau<2; resDau++){ | |
929 | Int_t indResDau=indFirstResDau+resDau; | |
930 | if(indResDau<0) return -1; | |
931 | TParticle* resdau=stack->Particle(indResDau); | |
932 | if(!resdau) return -1; | |
933 | Int_t pdgresdau=resdau->GetPdgCode(); | |
934 | if(TMath::Abs(pdgresdau)==321){ | |
935 | if(pdgD*pdgresdau>0) return -1; | |
936 | sumPxDau+=resdau->Px(); | |
937 | sumPyDau+=resdau->Py(); | |
938 | sumPzDau+=resdau->Pz(); | |
939 | nKaons++; | |
940 | arrayDauLab[nFoundKpi++]=indResDau; | |
941 | if(nFoundKpi>3) return -1; | |
942 | } | |
943 | if(TMath::Abs(pdgresdau)==211){ | |
944 | if(pdgD*pdgresdau<0) return -1; | |
945 | sumPxDau+=resdau->Px(); | |
946 | sumPyDau+=resdau->Py(); | |
947 | sumPzDau+=resdau->Pz(); | |
948 | nPions++; | |
949 | arrayDauLab[nFoundKpi++]=indResDau; | |
950 | if(nFoundKpi>3) return -1; | |
951 | } | |
952 | } | |
953 | }else{ | |
954 | return -1; | |
955 | } | |
956 | } | |
957 | if(nPions!=2) return -1; | |
958 | if(nKaons!=1) return -1; | |
959 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
960 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
961 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
962 | if(nDau==3) return 1; | |
963 | else if(nDau==2) return 2; | |
964 | } | |
965 | ||
966 | return -1; | |
967 | ||
968 | } | |
969 | ||
c8781624 | 970 | //____________________________________________________________________________ |
971 | Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){ | |
972 | // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases | |
973 | ||
974 | Int_t pdgD=mcPart->GetPdgCode(); | |
975 | if(TMath::Abs(pdgD)!=411) return -1; | |
976 | ||
977 | Int_t nDau=mcPart->GetNDaughters(); | |
978 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
979 | Int_t nKaons=0; | |
980 | Int_t nPions=0; | |
981 | Double_t sumPxDau=0.; | |
982 | Double_t sumPyDau=0.; | |
983 | Double_t sumPzDau=0.; | |
984 | Int_t nFoundKpi=0; | |
985 | ||
986 | if(nDau==3 || nDau==2){ | |
987 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
988 | Int_t indDau = labelFirstDau+iDau; | |
989 | if(indDau<0) return -1; | |
990 | AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau)); | |
991 | if(!dau) return -1; | |
992 | Int_t pdgdau=dau->GetPdgCode(); | |
993 | if(TMath::Abs(pdgdau)==321){ | |
994 | if(pdgD*pdgdau>0) return -1; | |
995 | nKaons++; | |
996 | sumPxDau+=dau->Px(); | |
997 | sumPyDau+=dau->Py(); | |
998 | sumPzDau+=dau->Pz(); | |
999 | arrayDauLab[nFoundKpi++]=indDau; | |
1000 | if(nFoundKpi>3) return -1; | |
1001 | }else if(TMath::Abs(pdgdau)==211){ | |
1002 | if(pdgD*pdgdau<0) return -1; | |
1003 | nPions++; | |
1004 | sumPxDau+=dau->Px(); | |
1005 | sumPyDau+=dau->Py(); | |
1006 | sumPzDau+=dau->Pz(); | |
1007 | arrayDauLab[nFoundKpi++]=indDau; | |
1008 | if(nFoundKpi>3) return -1; | |
1009 | }else if(TMath::Abs(pdgdau)==313){ | |
1010 | Int_t nResDau=dau->GetNDaughters(); | |
1011 | if(nResDau!=2) return -1; | |
1012 | Int_t indFirstResDau=dau->GetDaughter(0); | |
1013 | for(Int_t resDau=0; resDau<2; resDau++){ | |
1014 | Int_t indResDau=indFirstResDau+resDau; | |
1015 | if(indResDau<0) return -1; | |
1016 | AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau)); | |
1017 | if(!resdau) return -1; | |
1018 | Int_t pdgresdau=resdau->GetPdgCode(); | |
1019 | if(TMath::Abs(pdgresdau)==321){ | |
1020 | if(pdgD*pdgresdau>0) return -1; | |
1021 | sumPxDau+=resdau->Px(); | |
1022 | sumPyDau+=resdau->Py(); | |
1023 | sumPzDau+=resdau->Pz(); | |
1024 | nKaons++; | |
1025 | arrayDauLab[nFoundKpi++]=indResDau; | |
1026 | if(nFoundKpi>3) return -1; | |
1027 | } | |
1028 | if(TMath::Abs(pdgresdau)==211){ | |
1029 | if(pdgD*pdgresdau<0) return -1; | |
1030 | sumPxDau+=resdau->Px(); | |
1031 | sumPyDau+=resdau->Py(); | |
1032 | sumPzDau+=resdau->Pz(); | |
1033 | nPions++; | |
1034 | arrayDauLab[nFoundKpi++]=indResDau; | |
1035 | if(nFoundKpi>3) return -1; | |
1036 | } | |
1037 | } | |
1038 | }else{ | |
1039 | return -1; | |
1040 | } | |
1041 | } | |
1042 | if(nPions!=2) return -1; | |
1043 | if(nKaons!=1) return -1; | |
1044 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
1045 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
1046 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
1047 | if(nDau==3) return 1; | |
1048 | else if(nDau==2) return 2; | |
1049 | } | |
1050 | ||
1051 | return -1; | |
1052 | ||
5dc814ec | 1053 | } |
1054 | //____________________________________________________________________________ | |
1055 | Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ | |
1056 | // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case | |
1057 | ||
1058 | if(label<0) return -1; | |
1059 | TParticle* mcPart = stack->Particle(label); | |
1060 | Int_t pdgD=mcPart->GetPdgCode(); | |
1061 | if(TMath::Abs(pdgD)!=411) return -1; | |
1062 | ||
1063 | Int_t nDau=mcPart->GetNDaughters(); | |
1064 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
1065 | Int_t nKaons=0; | |
1066 | Int_t nPions=0; | |
1067 | Double_t sumPxDau=0.; | |
1068 | Double_t sumPyDau=0.; | |
1069 | Double_t sumPzDau=0.; | |
1070 | Int_t nFoundKpi=0; | |
1071 | Bool_t isPhi=kFALSE; | |
1072 | Bool_t isk0st=kFALSE; | |
1073 | ||
1074 | if(nDau==3 || nDau==2){ | |
1075 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
1076 | Int_t indDau = labelFirstDau+iDau; | |
1077 | if(indDau<0) return -1; | |
1078 | TParticle* dau=stack->Particle(indDau); | |
1079 | if(!dau) return -1; | |
1080 | Int_t pdgdau=dau->GetPdgCode(); | |
1081 | if(TMath::Abs(pdgdau)==321){ | |
1082 | nKaons++; | |
1083 | sumPxDau+=dau->Px(); | |
1084 | sumPyDau+=dau->Py(); | |
1085 | sumPzDau+=dau->Pz(); | |
1086 | arrayDauLab[nFoundKpi++]=indDau; | |
1087 | if(nFoundKpi>3) return -1; | |
1088 | }else if(TMath::Abs(pdgdau)==211){ | |
1089 | nPions++; | |
1090 | sumPxDau+=dau->Px(); | |
1091 | sumPyDau+=dau->Py(); | |
1092 | sumPzDau+=dau->Pz(); | |
1093 | arrayDauLab[nFoundKpi++]=indDau; | |
1094 | if(nFoundKpi>3) return -1; | |
1095 | }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){ | |
1096 | if(TMath::Abs(pdgdau)==313) isk0st=kTRUE; | |
1097 | if(TMath::Abs(pdgdau)==333) isPhi=kTRUE; | |
1098 | Int_t nResDau=dau->GetNDaughters(); | |
1099 | if(nResDau!=2) return -1; | |
1100 | Int_t indFirstResDau=dau->GetDaughter(0); | |
1101 | for(Int_t resDau=0; resDau<2; resDau++){ | |
1102 | Int_t indResDau=indFirstResDau+resDau; | |
1103 | if(indResDau<0) return -1; | |
1104 | TParticle* resdau=stack->Particle(indResDau); | |
1105 | if(!resdau) return -1; | |
1106 | Int_t pdgresdau=resdau->GetPdgCode(); | |
1107 | if(TMath::Abs(pdgresdau)==321){ | |
1108 | sumPxDau+=resdau->Px(); | |
1109 | sumPyDau+=resdau->Py(); | |
1110 | sumPzDau+=resdau->Pz(); | |
1111 | nKaons++; | |
1112 | arrayDauLab[nFoundKpi++]=indResDau; | |
1113 | if(nFoundKpi>3) return -1; | |
1114 | } | |
1115 | if(TMath::Abs(pdgresdau)==211){ | |
1116 | sumPxDau+=resdau->Px(); | |
1117 | sumPyDau+=resdau->Py(); | |
1118 | sumPzDau+=resdau->Pz(); | |
1119 | nPions++; | |
1120 | arrayDauLab[nFoundKpi++]=indResDau; | |
1121 | if(nFoundKpi>3) return -1; | |
1122 | } | |
1123 | } | |
1124 | }else{ | |
1125 | return -1; | |
1126 | } | |
1127 | } | |
1128 | if(nPions!=1) return -1; | |
1129 | if(nKaons!=2) return -1; | |
1130 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
1131 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
1132 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
1133 | if(nDau==3) return 3; | |
1134 | else if(nDau==2){ | |
1135 | if(isk0st) return 2; | |
1136 | if(isPhi) return 1; | |
1137 | } | |
1138 | } | |
1139 | ||
1140 | return -1; | |
1141 | ||
c8781624 | 1142 | } |
4d21d729 | 1143 | //____________________________________________________________________________ |
1144 | Int_t AliVertexingHFUtils::CheckDplusK0spiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ | |
1145 | // Checks the Dplus->V0+pion decay channel. Returns 1 if success, -1 otherwise | |
1146 | ||
1147 | if(label<0) return -1; | |
1148 | TParticle* mcPart = stack->Particle(label); | |
1149 | Int_t pdgD=mcPart->GetPdgCode(); | |
1150 | if(TMath::Abs(pdgD)!=411) return -1; | |
1151 | ||
1152 | Int_t nDau=mcPart->GetNDaughters(); | |
1153 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
1154 | Int_t nPions=0; | |
1155 | Double_t sumPxDau=0.; | |
1156 | Double_t sumPyDau=0.; | |
1157 | Double_t sumPzDau=0.; | |
1158 | Int_t nFoundpi=0; | |
1159 | ||
1160 | Int_t codeV0=-1; | |
1161 | if(nDau==2){ | |
1162 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
1163 | Int_t indDau = labelFirstDau+iDau; | |
1164 | if(indDau<0) return -1; | |
1165 | TParticle* dau=stack->Particle(indDau); | |
1166 | if(!dau) return -1; | |
1167 | Int_t pdgdau=dau->GetPdgCode(); | |
1168 | if(TMath::Abs(pdgdau)==211){ | |
1169 | nPions++; | |
1170 | sumPxDau+=dau->Px(); | |
1171 | sumPyDau+=dau->Py(); | |
1172 | sumPzDau+=dau->Pz(); | |
1173 | arrayDauLab[nFoundpi++]=indDau; | |
1174 | if(nFoundpi>3) return -1; | |
1175 | }else if(TMath::Abs(pdgdau)==311){ | |
1176 | codeV0=TMath::Abs(pdgdau); | |
1177 | TParticle* v0=dau; | |
1178 | if(codeV0==311){ | |
1179 | Int_t nK0Dau=dau->GetNDaughters(); | |
1180 | if(nK0Dau!=1) return -1; | |
1181 | Int_t indK0s=dau->GetDaughter(0); | |
1182 | if(indK0s<0) return -1; | |
1183 | v0=stack->Particle(indK0s); | |
1184 | if(!v0) return -1; | |
1185 | Int_t pdgK0sl=v0->GetPdgCode(); | |
1186 | codeV0=TMath::Abs(pdgK0sl); | |
1187 | } | |
1188 | Int_t nV0Dau=v0->GetNDaughters(); | |
1189 | if(nV0Dau!=2) return -1; | |
1190 | Int_t indFirstV0Dau=v0->GetDaughter(0); | |
1191 | for(Int_t v0Dau=0; v0Dau<2; v0Dau++){ | |
1192 | Int_t indV0Dau=indFirstV0Dau+v0Dau; | |
1193 | if(indV0Dau<0) return -1; | |
1194 | TParticle* v0dau=stack->Particle(indV0Dau); | |
1195 | if(!v0dau) return -1; | |
1196 | Int_t pdgv0dau=v0dau->GetPdgCode(); | |
1197 | if(TMath::Abs(pdgv0dau)==211){ | |
1198 | sumPxDau+=v0dau->Px(); | |
1199 | sumPyDau+=v0dau->Py(); | |
1200 | sumPzDau+=v0dau->Pz(); | |
1201 | nPions++; | |
1202 | arrayDauLab[nFoundpi++]=indV0Dau; | |
1203 | if(nFoundpi>3) return -1; | |
1204 | } | |
1205 | } | |
1206 | }else{ | |
1207 | return -1; | |
1208 | } | |
1209 | } | |
1210 | if(nPions!=3) return -1; | |
1211 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
1212 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
1213 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
1214 | if(codeV0==310) return 1; | |
1215 | } | |
1216 | return -1; | |
1217 | ||
1218 | } | |
1219 | ||
1220 | ||
ee9c9d9b | 1221 | //____________________________________________________________________________ |
1222 | Int_t AliVertexingHFUtils::CheckDsDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ | |
1223 | // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case | |
1224 | ||
1225 | if(label<0) return -1; | |
1226 | TParticle* mcPart = stack->Particle(label); | |
1227 | Int_t pdgD=mcPart->GetPdgCode(); | |
1228 | if(TMath::Abs(pdgD)!=431) return -1; | |
1229 | ||
1230 | Int_t nDau=mcPart->GetNDaughters(); | |
1231 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
1232 | Int_t nKaons=0; | |
1233 | Int_t nPions=0; | |
1234 | Double_t sumPxDau=0.; | |
1235 | Double_t sumPyDau=0.; | |
1236 | Double_t sumPzDau=0.; | |
1237 | Int_t nFoundKpi=0; | |
1238 | Bool_t isPhi=kFALSE; | |
1239 | Bool_t isk0st=kFALSE; | |
1240 | ||
1241 | if(nDau==3 || nDau==2){ | |
1242 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
1243 | Int_t indDau = labelFirstDau+iDau; | |
1244 | if(indDau<0) return -1; | |
1245 | TParticle* dau=stack->Particle(indDau); | |
1246 | if(!dau) return -1; | |
1247 | Int_t pdgdau=dau->GetPdgCode(); | |
1248 | if(TMath::Abs(pdgdau)==321){ | |
1249 | nKaons++; | |
1250 | sumPxDau+=dau->Px(); | |
1251 | sumPyDau+=dau->Py(); | |
1252 | sumPzDau+=dau->Pz(); | |
1253 | arrayDauLab[nFoundKpi++]=indDau; | |
1254 | if(nFoundKpi>3) return -1; | |
1255 | }else if(TMath::Abs(pdgdau)==211){ | |
1256 | nPions++; | |
1257 | sumPxDau+=dau->Px(); | |
1258 | sumPyDau+=dau->Py(); | |
1259 | sumPzDau+=dau->Pz(); | |
1260 | arrayDauLab[nFoundKpi++]=indDau; | |
1261 | if(nFoundKpi>3) return -1; | |
1262 | }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){ | |
1263 | if(TMath::Abs(pdgdau)==313) isk0st=kTRUE; | |
1264 | if(TMath::Abs(pdgdau)==333) isPhi=kTRUE; | |
1265 | Int_t nResDau=dau->GetNDaughters(); | |
1266 | if(nResDau!=2) return -1; | |
1267 | Int_t indFirstResDau=dau->GetDaughter(0); | |
1268 | for(Int_t resDau=0; resDau<2; resDau++){ | |
1269 | Int_t indResDau=indFirstResDau+resDau; | |
1270 | if(indResDau<0) return -1; | |
1271 | TParticle* resdau=stack->Particle(indResDau); | |
1272 | if(!resdau) return -1; | |
1273 | Int_t pdgresdau=resdau->GetPdgCode(); | |
1274 | if(TMath::Abs(pdgresdau)==321){ | |
1275 | sumPxDau+=resdau->Px(); | |
1276 | sumPyDau+=resdau->Py(); | |
1277 | sumPzDau+=resdau->Pz(); | |
1278 | nKaons++; | |
1279 | arrayDauLab[nFoundKpi++]=indResDau; | |
1280 | if(nFoundKpi>3) return -1; | |
1281 | } | |
1282 | if(TMath::Abs(pdgresdau)==211){ | |
1283 | sumPxDau+=resdau->Px(); | |
1284 | sumPyDau+=resdau->Py(); | |
1285 | sumPzDau+=resdau->Pz(); | |
1286 | nPions++; | |
1287 | arrayDauLab[nFoundKpi++]=indResDau; | |
1288 | if(nFoundKpi>3) return -1; | |
1289 | } | |
1290 | } | |
1291 | }else{ | |
1292 | return -1; | |
1293 | } | |
1294 | } | |
1295 | if(nPions!=1) return -1; | |
1296 | if(nKaons!=2) return -1; | |
1297 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
1298 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
1299 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
1300 | if(nDau==3) return 3; | |
1301 | else if(nDau==2){ | |
1302 | if(isk0st) return 2; | |
1303 | if(isPhi) return 1; | |
1304 | } | |
1305 | } | |
1306 | ||
1307 | return -1; | |
1308 | ||
1309 | } | |
4d21d729 | 1310 | //____________________________________________________________________________ |
1311 | Int_t AliVertexingHFUtils::CheckDsK0sKDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ | |
1312 | // Checks the Ds->K0s+S decay channel. Returns 1 in case of success, otherwise -1 | |
1313 | ||
1314 | if(label<0) return -1; | |
1315 | TParticle* mcPart = stack->Particle(label); | |
1316 | Int_t pdgD=mcPart->GetPdgCode(); | |
1317 | if(TMath::Abs(pdgD)!=431) return -1; | |
1318 | ||
1319 | Int_t nDau=mcPart->GetNDaughters(); | |
1320 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
1321 | Int_t nPions=0; | |
1322 | Int_t nKaons=0; | |
1323 | Double_t sumPxDau=0.; | |
1324 | Double_t sumPyDau=0.; | |
1325 | Double_t sumPzDau=0.; | |
1326 | Int_t nFoundKpi=0; | |
1327 | ||
1328 | Int_t codeV0=-1; | |
1329 | if(nDau==2){ | |
1330 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
1331 | Int_t indDau = labelFirstDau+iDau; | |
1332 | if(indDau<0) return -1; | |
1333 | TParticle* dau=stack->Particle(indDau); | |
1334 | if(!dau) return -1; | |
1335 | Int_t pdgdau=dau->GetPdgCode(); | |
1336 | if(TMath::Abs(pdgdau)==211){ | |
1337 | nPions++; | |
1338 | sumPxDau+=dau->Px(); | |
1339 | sumPyDau+=dau->Py(); | |
1340 | sumPzDau+=dau->Pz(); | |
1341 | arrayDauLab[nFoundKpi++]=indDau; | |
1342 | if(nFoundKpi>3) return -1; | |
1343 | }else if(TMath::Abs(pdgdau)==321){ | |
1344 | nKaons++; | |
1345 | sumPxDau+=dau->Px(); | |
1346 | sumPyDau+=dau->Py(); | |
1347 | sumPzDau+=dau->Pz(); | |
1348 | arrayDauLab[nFoundKpi++]=indDau; | |
1349 | if(nFoundKpi>3) return -1; | |
1350 | }else if(TMath::Abs(pdgdau)==311){ | |
1351 | codeV0=TMath::Abs(pdgdau); | |
1352 | TParticle* v0=dau; | |
1353 | if(codeV0==311){ | |
1354 | Int_t nK0Dau=dau->GetNDaughters(); | |
1355 | if(nK0Dau!=1) return -1; | |
1356 | Int_t indK0s=dau->GetDaughter(0); | |
1357 | if(indK0s<0) return -1; | |
1358 | v0=stack->Particle(indK0s); | |
1359 | if(!v0) return -1; | |
1360 | Int_t pdgK0sl=v0->GetPdgCode(); | |
1361 | codeV0=TMath::Abs(pdgK0sl); | |
1362 | } | |
1363 | Int_t nV0Dau=v0->GetNDaughters(); | |
1364 | if(nV0Dau!=2) return -1; | |
1365 | Int_t indFirstV0Dau=v0->GetDaughter(0); | |
1366 | for(Int_t v0Dau=0; v0Dau<2; v0Dau++){ | |
1367 | Int_t indV0Dau=indFirstV0Dau+v0Dau; | |
1368 | if(indV0Dau<0) return -1; | |
1369 | TParticle* v0dau=stack->Particle(indV0Dau); | |
1370 | if(!v0dau) return -1; | |
1371 | Int_t pdgv0dau=v0dau->GetPdgCode(); | |
1372 | if(TMath::Abs(pdgv0dau)==211){ | |
1373 | sumPxDau+=v0dau->Px(); | |
1374 | sumPyDau+=v0dau->Py(); | |
1375 | sumPzDau+=v0dau->Pz(); | |
1376 | nPions++; | |
1377 | arrayDauLab[nFoundKpi++]=indV0Dau; | |
1378 | if(nFoundKpi>3) return -1; | |
1379 | } | |
1380 | } | |
1381 | }else{ | |
1382 | return -1; | |
1383 | } | |
1384 | } | |
1385 | if(nPions!=2) return -1; | |
1386 | if(nKaons!=1) return -1; | |
1387 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
1388 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
1389 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
1390 | if(codeV0==310) return 1; | |
1391 | } | |
1392 | return -1; | |
1393 | ||
1394 | } | |
1395 | ||
1396 | ||
ee9c9d9b | 1397 | //____________________________________________________________________________ |
1398 | Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){ | |
1399 | // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case | |
1400 | ||
1401 | Int_t pdgD=mcPart->GetPdgCode(); | |
1402 | if(TMath::Abs(pdgD)!=431) return -1; | |
1403 | ||
1404 | Int_t nDau=mcPart->GetNDaughters(); | |
1405 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
1406 | Int_t nKaons=0; | |
1407 | Int_t nPions=0; | |
1408 | Double_t sumPxDau=0.; | |
1409 | Double_t sumPyDau=0.; | |
1410 | Double_t sumPzDau=0.; | |
1411 | Int_t nFoundKpi=0; | |
1412 | Bool_t isPhi=kFALSE; | |
1413 | Bool_t isk0st=kFALSE; | |
1414 | ||
1415 | if(nDau==3 || nDau==2){ | |
1416 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
1417 | Int_t indDau = labelFirstDau+iDau; | |
1418 | if(indDau<0) return -1; | |
1419 | AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau)); | |
1420 | if(!dau) return -1; | |
1421 | Int_t pdgdau=dau->GetPdgCode(); | |
1422 | if(TMath::Abs(pdgdau)==321){ | |
1423 | nKaons++; | |
1424 | sumPxDau+=dau->Px(); | |
1425 | sumPyDau+=dau->Py(); | |
1426 | sumPzDau+=dau->Pz(); | |
1427 | arrayDauLab[nFoundKpi++]=indDau; | |
1428 | if(nFoundKpi>3) return -1; | |
1429 | }else if(TMath::Abs(pdgdau)==211){ | |
1430 | nPions++; | |
1431 | sumPxDau+=dau->Px(); | |
1432 | sumPyDau+=dau->Py(); | |
1433 | sumPzDau+=dau->Pz(); | |
1434 | arrayDauLab[nFoundKpi++]=indDau; | |
1435 | if(nFoundKpi>3) return -1; | |
1436 | }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){ | |
1437 | if(TMath::Abs(pdgdau)==313) isk0st=kTRUE; | |
1438 | if(TMath::Abs(pdgdau)==333) isPhi=kTRUE; | |
1439 | Int_t nResDau=dau->GetNDaughters(); | |
1440 | if(nResDau!=2) return -1; | |
1441 | Int_t indFirstResDau=dau->GetDaughter(0); | |
1442 | for(Int_t resDau=0; resDau<2; resDau++){ | |
1443 | Int_t indResDau=indFirstResDau+resDau; | |
1444 | if(indResDau<0) return -1; | |
1445 | AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau)); | |
1446 | if(!resdau) return -1; | |
1447 | Int_t pdgresdau=resdau->GetPdgCode(); | |
1448 | if(TMath::Abs(pdgresdau)==321){ | |
1449 | sumPxDau+=resdau->Px(); | |
1450 | sumPyDau+=resdau->Py(); | |
1451 | sumPzDau+=resdau->Pz(); | |
1452 | nKaons++; | |
1453 | arrayDauLab[nFoundKpi++]=indResDau; | |
1454 | if(nFoundKpi>3) return -1; | |
1455 | } | |
1456 | if(TMath::Abs(pdgresdau)==211){ | |
1457 | sumPxDau+=resdau->Px(); | |
1458 | sumPyDau+=resdau->Py(); | |
1459 | sumPzDau+=resdau->Pz(); | |
1460 | nPions++; | |
1461 | arrayDauLab[nFoundKpi++]=indResDau; | |
1462 | if(nFoundKpi>3) return -1; | |
1463 | } | |
5dc814ec | 1464 | } |
ee9c9d9b | 1465 | }else{ |
1466 | return -1; | |
1467 | } | |
1468 | } | |
1469 | if(nPions!=1) return -1; | |
1470 | if(nKaons!=2) return -1; | |
1471 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
1472 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
1473 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
1474 | if(nDau==3) return 3; | |
1475 | else if(nDau==2){ | |
1476 | if(isk0st) return 2; | |
1477 | if(isPhi) return 1; | |
1478 | } | |
1479 | } | |
1480 | ||
1481 | return -1; | |
1482 | ||
1483 | } | |
1484 | //____________________________________________________________________________ | |
1485 | Int_t AliVertexingHFUtils::CheckDstarDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ | |
1486 | // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases | |
1487 | ||
1488 | if(label<0) return -1; | |
1489 | TParticle* mcPart = stack->Particle(label); | |
1490 | Int_t pdgD=mcPart->GetPdgCode(); | |
1491 | if(TMath::Abs(pdgD)!=413) return -1; | |
1492 | ||
1493 | Int_t nDau=mcPart->GetNDaughters(); | |
1494 | if(nDau!=2) return -1; | |
1495 | ||
1496 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
1497 | Int_t nKaons=0; | |
1498 | Int_t nPions=0; | |
1499 | Double_t sumPxDau=0.; | |
1500 | Double_t sumPyDau=0.; | |
1501 | Double_t sumPzDau=0.; | |
1502 | Int_t nFoundKpi=0; | |
1503 | ||
1504 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
1505 | Int_t indDau = labelFirstDau+iDau; | |
1506 | if(indDau<0) return -1; | |
1507 | TParticle* dau=stack->Particle(indDau); | |
1508 | if(!dau) return -1; | |
1509 | Int_t pdgdau=dau->GetPdgCode(); | |
1510 | if(TMath::Abs(pdgdau)==421){ | |
1511 | Int_t nResDau=dau->GetNDaughters(); | |
1512 | if(nResDau!=2) return -1; | |
1513 | Int_t indFirstResDau=dau->GetDaughter(0); | |
1514 | for(Int_t resDau=0; resDau<2; resDau++){ | |
1515 | Int_t indResDau=indFirstResDau+resDau; | |
1516 | if(indResDau<0) return -1; | |
1517 | TParticle* resdau=stack->Particle(indResDau); | |
1518 | if(!resdau) return -1; | |
1519 | Int_t pdgresdau=resdau->GetPdgCode(); | |
1520 | if(TMath::Abs(pdgresdau)==321){ | |
1521 | if(pdgD*pdgresdau>0) return -1; | |
1522 | sumPxDau+=resdau->Px(); | |
1523 | sumPyDau+=resdau->Py(); | |
1524 | sumPzDau+=resdau->Pz(); | |
1525 | nKaons++; | |
1526 | arrayDauLab[nFoundKpi++]=indResDau; | |
1527 | if(nFoundKpi>3) return -1; | |
1528 | } | |
1529 | if(TMath::Abs(pdgresdau)==211){ | |
1530 | if(pdgD*pdgresdau<0) return -1; | |
1531 | sumPxDau+=resdau->Px(); | |
1532 | sumPyDau+=resdau->Py(); | |
1533 | sumPzDau+=resdau->Pz(); | |
1534 | nPions++; | |
1535 | arrayDauLab[nFoundKpi++]=indResDau; | |
1536 | if(nFoundKpi>3) return -1; | |
1537 | } | |
1538 | } | |
1539 | }else if(TMath::Abs(pdgdau)==211){ | |
1540 | if(pdgD*pdgdau<0) return -1; | |
1541 | nPions++; | |
1542 | sumPxDau+=dau->Px(); | |
1543 | sumPyDau+=dau->Py(); | |
1544 | sumPzDau+=dau->Pz(); | |
1545 | arrayDauLab[nFoundKpi++]=indDau; | |
1546 | if(nFoundKpi>3) return -1; | |
1547 | } | |
1548 | } | |
1549 | ||
1550 | if(nPions!=2) return -1; | |
1551 | if(nKaons!=1) return -1; | |
1552 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
1553 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
1554 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
1555 | return 1; | |
1556 | ||
1557 | } | |
1558 | ||
1559 | //____________________________________________________________________________ | |
1560 | Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){ | |
1561 | // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases | |
1562 | ||
1563 | Int_t pdgD=mcPart->GetPdgCode(); | |
1564 | if(TMath::Abs(pdgD)!=413) return -1; | |
1565 | ||
1566 | Int_t nDau=mcPart->GetNDaughters(); | |
1567 | if(nDau!=2) return -1; | |
1568 | ||
1569 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
1570 | Int_t nKaons=0; | |
1571 | Int_t nPions=0; | |
1572 | Double_t sumPxDau=0.; | |
1573 | Double_t sumPyDau=0.; | |
1574 | Double_t sumPzDau=0.; | |
1575 | Int_t nFoundKpi=0; | |
1576 | ||
1577 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
1578 | Int_t indDau = labelFirstDau+iDau; | |
1579 | if(indDau<0) return -1; | |
1580 | AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau)); | |
1581 | if(!dau) return -1; | |
1582 | Int_t pdgdau=dau->GetPdgCode(); | |
1583 | if(TMath::Abs(pdgdau)==421){ | |
1584 | Int_t nResDau=dau->GetNDaughters(); | |
1585 | if(nResDau!=2) return -1; | |
1586 | Int_t indFirstResDau=dau->GetDaughter(0); | |
1587 | for(Int_t resDau=0; resDau<2; resDau++){ | |
1588 | Int_t indResDau=indFirstResDau+resDau; | |
1589 | if(indResDau<0) return -1; | |
1590 | AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau)); | |
1591 | if(!resdau) return -1; | |
1592 | Int_t pdgresdau=resdau->GetPdgCode(); | |
1593 | if(TMath::Abs(pdgresdau)==321){ | |
1594 | if(pdgD*pdgresdau>0) return -1; | |
1595 | sumPxDau+=resdau->Px(); | |
1596 | sumPyDau+=resdau->Py(); | |
1597 | sumPzDau+=resdau->Pz(); | |
1598 | nKaons++; | |
1599 | arrayDauLab[nFoundKpi++]=indResDau; | |
1600 | if(nFoundKpi>3) return -1; | |
1601 | } | |
1602 | if(TMath::Abs(pdgresdau)==211){ | |
1603 | if(pdgD*pdgresdau<0) return -1; | |
1604 | sumPxDau+=resdau->Px(); | |
1605 | sumPyDau+=resdau->Py(); | |
1606 | sumPzDau+=resdau->Pz(); | |
1607 | nPions++; | |
1608 | arrayDauLab[nFoundKpi++]=indResDau; | |
1609 | if(nFoundKpi>3) return -1; | |
1610 | } | |
1611 | } | |
1612 | }else if(TMath::Abs(pdgdau)==211){ | |
1613 | if(pdgD*pdgdau<0) return -1; | |
1614 | nPions++; | |
1615 | sumPxDau+=dau->Px(); | |
1616 | sumPyDau+=dau->Py(); | |
1617 | sumPzDau+=dau->Pz(); | |
1618 | arrayDauLab[nFoundKpi++]=indDau; | |
1619 | if(nFoundKpi>3) return -1; | |
1620 | } | |
1621 | } | |
1622 | ||
1623 | if(nPions!=2) return -1; | |
1624 | if(nKaons!=1) return -1; | |
1625 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
1626 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
1627 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
1628 | return 1; | |
1629 | ||
1630 | } | |
1631 | ||
5dc814ec | 1632 | //____________________________________________________________________________ |
1633 | Int_t AliVertexingHFUtils::CheckLcpKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ | |
1634 | // Checks the Lc->pKpi decay channel. Returns 1 for non-resonant decays and 2, 3 or 4 for resonant ones, -1 in other cases | |
1635 | ||
1636 | if(label<0) return -1; | |
1637 | TParticle* mcPart = stack->Particle(label); | |
1638 | Int_t pdgD=mcPart->GetPdgCode(); | |
1639 | if(TMath::Abs(pdgD)!=4122) return -1; | |
1640 | ||
1641 | Int_t nDau=mcPart->GetNDaughters(); | |
1642 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
1643 | Int_t nKaons=0; | |
1644 | Int_t nPions=0; | |
1645 | Int_t nProtons=0; | |
1646 | Double_t sumPxDau=0.; | |
1647 | Double_t sumPyDau=0.; | |
1648 | Double_t sumPzDau=0.; | |
1649 | Int_t nFoundpKpi=0; | |
1650 | ||
1651 | Int_t codeRes=-1; | |
1652 | if(nDau==3 || nDau==2){ | |
1653 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
1654 | Int_t indDau = labelFirstDau+iDau; | |
1655 | if(indDau<0) return -1; | |
1656 | TParticle* dau=stack->Particle(indDau); | |
1657 | if(!dau) return -1; | |
1658 | Int_t pdgdau=dau->GetPdgCode(); | |
1659 | if(TMath::Abs(pdgdau)==321){ | |
1660 | nKaons++; | |
1661 | sumPxDau+=dau->Px(); | |
1662 | sumPyDau+=dau->Py(); | |
1663 | sumPzDau+=dau->Pz(); | |
1664 | arrayDauLab[nFoundpKpi++]=indDau; | |
1665 | if(nFoundpKpi>3) return -1; | |
1666 | }else if(TMath::Abs(pdgdau)==211){ | |
1667 | nPions++; | |
1668 | sumPxDau+=dau->Px(); | |
1669 | sumPyDau+=dau->Py(); | |
1670 | sumPzDau+=dau->Pz(); | |
1671 | arrayDauLab[nFoundpKpi++]=indDau; | |
1672 | if(nFoundpKpi>3) return -1; | |
1673 | }else if(TMath::Abs(pdgdau)==2212){ | |
1674 | nProtons++; | |
1675 | sumPxDau+=dau->Px(); | |
1676 | sumPyDau+=dau->Py(); | |
1677 | sumPzDau+=dau->Pz(); | |
1678 | arrayDauLab[nFoundpKpi++]=indDau; | |
1679 | if(nFoundpKpi>3) return -1; | |
1680 | }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 || | |
1681 | TMath::Abs(pdgdau)==2224){ | |
1682 | codeRes=TMath::Abs(pdgdau); | |
1683 | Int_t nResDau=dau->GetNDaughters(); | |
1684 | if(nResDau!=2) return -1; | |
1685 | Int_t indFirstResDau=dau->GetDaughter(0); | |
1686 | for(Int_t resDau=0; resDau<2; resDau++){ | |
1687 | Int_t indResDau=indFirstResDau+resDau; | |
1688 | if(indResDau<0) return -1; | |
1689 | TParticle* resdau=stack->Particle(indResDau); | |
1690 | if(!resdau) return -1; | |
1691 | Int_t pdgresdau=resdau->GetPdgCode(); | |
1692 | if(TMath::Abs(pdgresdau)==321){ | |
1693 | sumPxDau+=resdau->Px(); | |
1694 | sumPyDau+=resdau->Py(); | |
1695 | sumPzDau+=resdau->Pz(); | |
1696 | nKaons++; | |
1697 | arrayDauLab[nFoundpKpi++]=indResDau; | |
1698 | if(nFoundpKpi>3) return -1; | |
1699 | }else if(TMath::Abs(pdgresdau)==211){ | |
1700 | sumPxDau+=resdau->Px(); | |
1701 | sumPyDau+=resdau->Py(); | |
1702 | sumPzDau+=resdau->Pz(); | |
1703 | nPions++; | |
1704 | arrayDauLab[nFoundpKpi++]=indResDau; | |
1705 | if(nFoundpKpi>3) return -1; | |
1706 | }else if(TMath::Abs(pdgresdau)==2212){ | |
1707 | sumPxDau+=resdau->Px(); | |
1708 | sumPyDau+=resdau->Py(); | |
1709 | sumPzDau+=resdau->Pz(); | |
1710 | nProtons++; | |
1711 | arrayDauLab[nFoundpKpi++]=indResDau; | |
1712 | if(nFoundpKpi>3) return -1; | |
1713 | } | |
1714 | } | |
1715 | }else{ | |
1716 | return -1; | |
1717 | } | |
1718 | } | |
1719 | if(nPions!=1) return -1; | |
1720 | if(nKaons!=1) return -1; | |
1721 | if(nProtons!=1) return -1; | |
1722 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
1723 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
1724 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
1725 | if(nDau==3) return 1; | |
1726 | else if(nDau==2){ | |
1727 | if(codeRes==313) return 2; | |
1728 | else if(codeRes==2224) return 3; | |
1729 | else if(codeRes==3124) return 4; | |
1730 | } | |
1731 | } | |
1732 | return -1; | |
1733 | ||
1734 | } | |
1735 | ||
1736 | ||
1737 | //____________________________________________________________________________ | |
1738 | Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ | |
1739 | // Checks the Lc->V0+bachelor decay channel. Returns 1 for pK0s, 2 for piLambda, 3 for pK0l -1 in other cases | |
1740 | ||
1741 | if(label<0) return -1; | |
1742 | TParticle* mcPart = stack->Particle(label); | |
1743 | Int_t pdgD=mcPart->GetPdgCode(); | |
1744 | if(TMath::Abs(pdgD)!=4122) return -1; | |
1745 | ||
1746 | Int_t nDau=mcPart->GetNDaughters(); | |
1747 | Int_t labelFirstDau = mcPart->GetDaughter(0); | |
1748 | Int_t nPions=0; | |
1749 | Int_t nProtons=0; | |
1750 | Double_t sumPxDau=0.; | |
1751 | Double_t sumPyDau=0.; | |
1752 | Double_t sumPzDau=0.; | |
1753 | Int_t nFoundppi=0; | |
1754 | ||
1755 | Int_t codeV0=-1; | |
1756 | if(nDau==2){ | |
1757 | for(Int_t iDau=0; iDau<nDau; iDau++){ | |
1758 | Int_t indDau = labelFirstDau+iDau; | |
1759 | if(indDau<0) return -1; | |
1760 | TParticle* dau=stack->Particle(indDau); | |
1761 | if(!dau) return -1; | |
1762 | Int_t pdgdau=dau->GetPdgCode(); | |
1763 | if(TMath::Abs(pdgdau)==211){ | |
1764 | nPions++; | |
1765 | sumPxDau+=dau->Px(); | |
1766 | sumPyDau+=dau->Py(); | |
1767 | sumPzDau+=dau->Pz(); | |
1768 | arrayDauLab[nFoundppi++]=indDau; | |
1769 | if(nFoundppi>3) return -1; | |
1770 | }else if(TMath::Abs(pdgdau)==2212){ | |
1771 | nProtons++; | |
1772 | sumPxDau+=dau->Px(); | |
1773 | sumPyDau+=dau->Py(); | |
1774 | sumPzDau+=dau->Pz(); | |
1775 | arrayDauLab[nFoundppi++]=indDau; | |
1776 | if(nFoundppi>3) return -1; | |
1777 | }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){ | |
1778 | codeV0=TMath::Abs(pdgdau); | |
1779 | TParticle* v0=dau; | |
1780 | if(codeV0==311){ | |
1781 | Int_t nK0Dau=dau->GetNDaughters(); | |
1782 | if(nK0Dau!=1) return -1; | |
1783 | Int_t indK0s=dau->GetDaughter(0); | |
1784 | if(indK0s<0) return -1; | |
1785 | v0=stack->Particle(indK0s); | |
1786 | if(!v0) return -1; | |
1787 | Int_t pdgK0sl=v0->GetPdgCode(); | |
1788 | codeV0=TMath::Abs(pdgK0sl); | |
1789 | } | |
1790 | Int_t nV0Dau=v0->GetNDaughters(); | |
1791 | if(nV0Dau!=2) return -1; | |
1792 | Int_t indFirstV0Dau=v0->GetDaughter(0); | |
1793 | for(Int_t v0Dau=0; v0Dau<2; v0Dau++){ | |
1794 | Int_t indV0Dau=indFirstV0Dau+v0Dau; | |
1795 | if(indV0Dau<0) return -1; | |
1796 | TParticle* v0dau=stack->Particle(indV0Dau); | |
1797 | if(!v0dau) return -1; | |
1798 | Int_t pdgv0dau=v0dau->GetPdgCode(); | |
1799 | if(TMath::Abs(pdgv0dau)==211){ | |
1800 | sumPxDau+=v0dau->Px(); | |
1801 | sumPyDau+=v0dau->Py(); | |
1802 | sumPzDau+=v0dau->Pz(); | |
1803 | nPions++; | |
1804 | arrayDauLab[nFoundppi++]=indV0Dau; | |
1805 | if(nFoundppi>3) return -1; | |
1806 | }else if(TMath::Abs(pdgv0dau)==2212){ | |
1807 | sumPxDau+=v0dau->Px(); | |
1808 | sumPyDau+=v0dau->Py(); | |
1809 | sumPzDau+=v0dau->Pz(); | |
1810 | nProtons++; | |
1811 | arrayDauLab[nFoundppi++]=indV0Dau; | |
1812 | if(nFoundppi>3) return -1; | |
1813 | } | |
1814 | } | |
1815 | }else{ | |
1816 | return -1; | |
1817 | } | |
1818 | } | |
1819 | if(nPions!=2) return -1; | |
1820 | if(nProtons!=1) return -1; | |
1821 | if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; | |
1822 | if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; | |
1823 | if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; | |
1824 | if(codeV0==310) return 1; | |
1825 | else if(codeV0==3122) return 2; | |
1826 | } | |
1827 | return -1; | |
1828 | ||
1829 | } | |
1830 | ||
1831 |