]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliVertexingHFUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' into TPCdev
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliVertexingHFUtils.cxx
CommitLineData
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
42ClassImp(AliVertexingHFUtils)
43
44//______________________________________________________________________
45AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
46 fK(1),
47 fSubRes(1.),
48 fMinEtaForTracklets(-1.),
49 fMaxEtaForTracklets(1.)
50{
51 // Default contructor
52}
53
54
55//______________________________________________________________________
56AliVertexingHFUtils::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//______________________________________________________________________
68void 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//______________________________________________________________________
86Double_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//______________________________________________________________________
100Double_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//______________________________________________________________________
106Double_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//______________________________________________________________________
149Double_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//______________________________________________________________________
162Double_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//______________________________________________________________________
169Double_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//______________________________________________________________________
177Double_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//______________________________________________________________________
185Int_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 198Int_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//______________________________________________________________________
211Int_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//______________________________________________________________________
226Int_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 241void 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//____________________________________________________________________________
318Double_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//____________________________________________________________________________
361Double_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//____________________________________________________________________________
447Double_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//______________________________________________________________________
472TString 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//_____________________________________________________________________
489void 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//----------------------------------------------------------------------
521Bool_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//____________________________________________________________________________
531Bool_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 542Bool_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 562Int_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//____________________________________________________________________________
594Int_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//____________________________________________________________________________
702Int_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}