Added protection for null pointer
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliVertexingHFUtils.cxx
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>
17 #include <TRandom.h>
18 #include <TProfile.h>
19 #include <TClonesArray.h>
20 #include <TH1F.h>
21 #include <TH2F.h>
22 #include <TF1.h>
23 #include <TParticle.h>
24 #include "AliStack.h"
25 #include "AliAODEvent.h"
26 #include "AliAODMCHeader.h"
27 #include "AliGenEventHeader.h"
28 #include "AliAODMCParticle.h"
29 #include "AliAODRecoDecayHF.h"
30 #include "AliVertexingHFUtils.h"
31
32 /* $Id$ */
33
34 ///////////////////////////////////////////////////////////////////
35 //                                                               //
36 // Class with functions useful for different D2H analyses        //
37 // - event plane resolution                                      //
38 // - <pt> calculation with side band subtraction                 //
39 // - tracklet multiplicity calculation                            //
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
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 }
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 }
199 //______________________________________________________________________
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 }
242 //______________________________________________________________________
243 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){
244
245   // Compute <pt> from 2D histogram M vs pt
246
247   //Make 2D histos in the desired pt range
248   Int_t start=hMassD->FindBin(ptmin);
249   Int_t end=hMassD->FindBin(ptmax)-1;
250   const Int_t nx=end-start;
251   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()));
252   for(Int_t ix=start;ix<end;ix++){
253     for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
254       hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
255       hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
256     }
257   }
258
259   Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
260   Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
261   Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
262   Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
263   Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
264   Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
265   //  printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
266
267   Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
268   Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2);
269   Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
270   Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
271   Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
272   Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
273   Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
274   Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1);
275   Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
276   Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
277   //  printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
278
279   Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
280   Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
281   Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
282   //  printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
283
284   TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
285   TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
286   TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
287
288   hMptBkgLo->Rebin(rebin);
289   hMptBkgHi->Rebin(rebin);
290   hMptSigReg->Rebin(rebin);
291
292   hMptBkgLo->Sumw2();
293   hMptBkgHi->Sumw2();
294   TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
295   hMptBkgLoScal->Scale(bkgSig/bkgLow);
296   TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
297   hMptBkgHiScal->Scale(bkgSig/bkgHi);
298
299   TH1F* hMptBkgAver=0x0;
300   hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
301   hMptBkgAver->Add(hMptBkgHiScal);
302   hMptBkgAver->Scale(0.5);
303   TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
304   hMptSig->Add(hMptBkgAver,-1.);   
305  
306   averagePt = hMptSig->GetMean();
307   errorPt = hMptSig->GetMeanError();
308
309   delete hMptBkgLo;
310   delete hMptBkgHi;
311   delete hMptSigReg;
312   delete hMptBkgLoScal;
313   delete hMptBkgHiScal;
314   delete hMptBkgAver;
315   delete hMassDpt;
316   delete hMptSig;
317
318 }
319 //____________________________________________________________________________
320 Bool_t AliVertexingHFUtils::CheckT0TriggerFired(AliAODEvent* aodEv){
321   // check if T0VTX trigger was fired, based on a workaround suggested by Alla
322   const Double32_t *mean = aodEv->GetT0TOF();
323   if(mean && mean[0]<9999.) return kTRUE;
324   else return kFALSE;
325 }
326 //____________________________________________________________________________
327 Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
328   // true impact parameter calculation for Dzero
329
330   if(!partD || !arrayMC || !mcHeader) return 99999.;
331   Int_t code=partD->GetPdgCode();
332   if(TMath::Abs(code)!=421) return 99999.;
333
334   Double_t vtxTrue[3];
335   mcHeader->GetVertex(vtxTrue);
336   Double_t origD[3];
337   partD->XvYvZv(origD);
338   Short_t charge=partD->Charge();
339   Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
340   for(Int_t iDau=0; iDau<2; iDau++){
341     pXdauTrue[iDau]=0.;
342     pYdauTrue[iDau]=0.;
343     pZdauTrue[iDau]=0.;
344   }
345
346   Int_t nDau=partD->GetNDaughters();
347   Int_t labelFirstDau = partD->GetDaughter(0); 
348   if(nDau==2){
349     for(Int_t iDau=0; iDau<2; iDau++){
350       Int_t ind = labelFirstDau+iDau;
351       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
352       if(!part){
353         printf("Daughter particle not found in MC array");
354         return 99999.;
355       } 
356       pXdauTrue[iDau]=part->Px();
357       pYdauTrue[iDau]=part->Py();
358       pZdauTrue[iDau]=part->Pz();
359     }
360   }else{
361     return 99999.;
362   }
363
364   Double_t d0dummy[2]={0.,0.};
365   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
366   return aodDvsMC.ImpParXY();
367
368 }
369 //____________________________________________________________________________
370 Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
371   // true impact parameter calculation for Dplus
372
373   if(!partD || !arrayMC || !mcHeader) return 99999.;
374   Int_t code=partD->GetPdgCode();
375   if(TMath::Abs(code)!=411) return 99999.;
376
377   Double_t vtxTrue[3];
378   mcHeader->GetVertex(vtxTrue);
379   Double_t origD[3];
380   partD->XvYvZv(origD);
381   Short_t charge=partD->Charge();
382   Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
383   for(Int_t iDau=0; iDau<3; iDau++){
384     pXdauTrue[iDau]=0.;
385     pYdauTrue[iDau]=0.;
386     pZdauTrue[iDau]=0.;
387   }
388
389   Int_t nDau=partD->GetNDaughters();
390   Int_t labelFirstDau = partD->GetDaughter(0); 
391   if(nDau==3){
392     for(Int_t iDau=0; iDau<3; iDau++){
393       Int_t ind = labelFirstDau+iDau;
394       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
395       if(!part){
396         printf("Daughter particle not found in MC array");
397         return 99999.;
398       } 
399       pXdauTrue[iDau]=part->Px();
400       pYdauTrue[iDau]=part->Py();
401       pZdauTrue[iDau]=part->Pz();
402     }
403   }else if(nDau==2){
404     Int_t theDau=0;
405     for(Int_t iDau=0; iDau<2; iDau++){
406       Int_t ind = labelFirstDau+iDau;
407       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
408       if(!part){
409         printf("Daughter particle not found in MC array");
410         return 99999.;
411       } 
412       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
413       if(pdgCode==211 || pdgCode==321){
414         pXdauTrue[theDau]=part->Px();
415         pYdauTrue[theDau]=part->Py();
416         pZdauTrue[theDau]=part->Pz();
417         ++theDau;
418       }else{
419         Int_t nDauRes=part->GetNDaughters();
420         if(nDauRes==2){
421           Int_t labelFirstDauRes = part->GetDaughter(0);        
422           for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
423             Int_t indDR = labelFirstDauRes+iDauRes;
424             AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
425             if(!partDR){
426               printf("Daughter particle not found in MC array");
427               return 99999.;
428             } 
429             
430             Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
431             if(pdgCodeDR==211 || pdgCodeDR==321){
432               pXdauTrue[theDau]=partDR->Px();
433               pYdauTrue[theDau]=partDR->Py();
434               pZdauTrue[theDau]=partDR->Pz();
435               ++theDau;
436             }
437           }
438         }
439       }
440     }
441     if(theDau!=3){
442       printf("Wrong number of decay prongs");
443       return 99999.;
444     }
445   }
446
447   Double_t d0dummy[3]={0.,0.,0.};
448   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
449   return aodDvsMC.ImpParXY();
450
451 }
452
453
454
455 //____________________________________________________________________________
456 Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
457   //
458   // Correct the number of accepted tracklets based on a TProfile Hist
459   //
460   //
461
462   if(TMath::Abs(vtxZ)>10.0){
463     //    printf("ERROR: Z vertex out of range for correction of multiplicity\n");
464     return uncorrectedNacc;
465   }
466
467   if(!estimatorAvg){
468     printf("ERROR: Missing TProfile for correction of multiplicity\n");
469     return uncorrectedNacc;
470   }
471
472   Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
473
474   Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
475
476   Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
477
478   return correctedNacc;
479 }
480 //______________________________________________________________________
481 TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
482   // get the name of the generator that produced a given particle
483
484   Int_t nsumpart=0;
485   TList *lh=header->GetCocktailHeaders();
486   Int_t nh=lh->GetEntries();
487   for(Int_t i=0;i<nh;i++){
488     AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
489     TString genname=gh->GetName();
490     Int_t npart=gh->NProduced();
491     if(label>=nsumpart && label<(nsumpart+npart)) return genname;
492     nsumpart+=npart;
493   }
494   TString empty="";
495   return empty;
496 }
497 //_____________________________________________________________________
498 void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
499
500   // method to check if a track comes from a given generator
501
502   Int_t lab=TMath::Abs(track->GetLabel());
503   nameGen=GetGenerator(lab,header);
504   
505   //  Int_t countControl=0;
506   
507   while(nameGen.IsWhitespace()){
508     AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
509     if(!mcpart){
510       printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
511       break;
512     }
513     Int_t mother = mcpart->GetMother();
514     if(mother<0){
515       printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
516       break;
517     }
518     lab=mother;
519     nameGen=GetGenerator(mother,header);
520     // countControl++;
521     // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
522     //   printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
523     //   break;
524     // }
525   }
526   
527   return;
528 }
529 //----------------------------------------------------------------------
530 Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
531   // method to check if a track comes from the signal event or from the underlying Hijing event
532   TString nameGen;
533   GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
534   
535   if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
536   
537   return kTRUE;
538 }
539 //____________________________________________________________________________
540 Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
541   // method to check if a D meson candidate comes from the signal event or from the underlying Hijing event
542
543   Int_t nprongs=cand->GetNProngs();
544   for(Int_t i=0;i<nprongs;i++){
545     AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i);
546     if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
547   }
548   return kFALSE;
549 }
550 //____________________________________________________________________________
551 Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
552   // method to check if a cascade candidate comes from the signal event or from the underlying Hijing event
553
554   AliAODTrack* bach = cand->GetBachelor();
555   if(IsTrackInjected(bach, header, arrayMC)) {
556     AliDebug(2, "Bachelor is injected, the whole candidate is then injected");
557     return kTRUE;
558   }
559   AliAODv0* v0 = cand->Getv0();
560   Int_t nprongs = v0->GetNProngs();
561   for(Int_t i = 0; i < nprongs; i++){
562     AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i);
563     if(IsTrackInjected(daugh,header,arrayMC)) {
564       AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i));
565       return kTRUE;
566     }
567   }
568   return kFALSE;
569 }
570 //____________________________________________________________________________
571 Int_t AliVertexingHFUtils::CheckOrigin(AliStack* stack, TParticle *mcPart, Bool_t searchUpToQuark){
572   // checking whether the mother of the particles come from a charm or a bottom quark
573
574   Int_t pdgGranma = 0;
575   Int_t mother = 0;
576   mother = mcPart->GetFirstMother();
577   Int_t istep = 0;
578   Int_t abspdgGranma =0;
579   Bool_t isFromB=kFALSE;
580   Bool_t isQuarkFound=kFALSE;
581   while (mother >0 ){
582     istep++;
583     TParticle* mcGranma = stack->Particle(mother);
584     if (mcGranma){
585       pdgGranma = mcGranma->GetPdgCode();
586       abspdgGranma = TMath::Abs(pdgGranma);
587       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
588         isFromB=kTRUE;
589       }
590       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
591       mother = mcGranma->GetFirstMother();
592     }else{
593       printf("CheckOrigin: Failed casting the mother particle!");
594       break;
595     }
596   }
597   if(searchUpToQuark && !isQuarkFound) return 0;
598   if(isFromB) return 5;
599   else return 4;
600
601 }
602 //____________________________________________________________________________
603 Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
604   // checking whether the mother of the particles come from a charm or a bottom quark
605
606   Int_t pdgGranma = 0;
607   Int_t mother = 0;
608   mother = mcPart->GetMother();
609   Int_t istep = 0;
610   Int_t abspdgGranma =0;
611   Bool_t isFromB=kFALSE;
612   Bool_t isQuarkFound=kFALSE;
613   while (mother >0 ){
614     istep++;
615     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
616     if (mcGranma){
617       pdgGranma = mcGranma->GetPdgCode();
618       abspdgGranma = TMath::Abs(pdgGranma);
619       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
620         isFromB=kTRUE;
621       }
622       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
623       mother = mcGranma->GetMother();
624     }else{
625       printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
626       break;
627     }
628   }
629   if(searchUpToQuark && !isQuarkFound) return 0;
630   if(isFromB) return 5;
631   else return 4;
632
633 }
634
635 //____________________________________________________________________________
636 Int_t AliVertexingHFUtils::CheckD0Decay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
637   // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
638
639   if(label<0) return -1;
640   TParticle* mcPart = stack->Particle(label);
641   Int_t pdgD=mcPart->GetPdgCode();
642   if(TMath::Abs(pdgD)!=421) return -1;
643   
644   Int_t nDau=mcPart->GetNDaughters();
645   
646   if(nDau==2){
647     Int_t daughter0 = mcPart->GetDaughter(0);
648     Int_t daughter1 = mcPart->GetDaughter(1);
649     TParticle* mcPartDaughter0 = stack->Particle(daughter0);
650     TParticle* mcPartDaughter1 = stack->Particle(daughter1);
651     if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
652     arrayDauLab[0]=daughter0;
653     arrayDauLab[1]=daughter1;
654     Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
655     Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
656     if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
657        !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
658       return -1;
659     }
660     if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
661     if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
662     if((pdgCode0*pdgCode1)>0) return -1;
663     Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
664     Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
665     Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
666     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
667     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
668     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
669     return 1;
670   }
671   
672   if(nDau==3 || nDau==4){
673     Int_t nKaons=0;
674     Int_t nPions=0;
675     Double_t sumPxDau=0.;
676     Double_t sumPyDau=0.;
677     Double_t sumPzDau=0.;
678     Int_t nFoundKpi=0;
679     Int_t labelFirstDau = mcPart->GetDaughter(0);
680     for(Int_t iDau=0; iDau<nDau; iDau++){
681       Int_t indDau = labelFirstDau+iDau;
682       if(indDau<0) return -1;
683       TParticle* dau=stack->Particle(indDau);
684       if(!dau) return -1;
685       Int_t pdgdau=dau->GetPdgCode();
686       if(TMath::Abs(pdgdau)==321){
687         if(pdgD>0 && pdgdau>0) return -1;
688         if(pdgD<0 && pdgdau<0) return -1;
689         nKaons++;
690         sumPxDau+=dau->Px();
691         sumPyDau+=dau->Py();
692         sumPzDau+=dau->Pz();
693         arrayDauLab[nFoundKpi++]=indDau;
694         if(nFoundKpi>4) return -1;
695       }else if(TMath::Abs(pdgdau)==211){
696         nPions++;
697         sumPxDau+=dau->Px();
698         sumPyDau+=dau->Py();
699         sumPzDau+=dau->Pz();
700         arrayDauLab[nFoundKpi++]=indDau;
701         if(nFoundKpi>4) return -1;
702       }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
703         Int_t nResDau=dau->GetNDaughters();
704         if(nResDau!=2) return -1;
705         Int_t indFirstResDau=dau->GetDaughter(0);
706         for(Int_t resDau=0; resDau<2; resDau++){
707           Int_t indResDau=indFirstResDau+resDau;
708           if(indResDau<0) return -1;
709           TParticle* resdau=stack->Particle(indResDau);
710           if(!resdau) return -1;
711           Int_t pdgresdau=resdau->GetPdgCode();
712           if(TMath::Abs(pdgresdau)==321){
713             if(pdgD>0 && pdgresdau>0) return -1;
714             if(pdgD<0 && pdgresdau<0) return -1;
715             nKaons++;
716             sumPxDau+=resdau->Px();
717             sumPyDau+=resdau->Py();
718             sumPzDau+=resdau->Pz();
719             arrayDauLab[nFoundKpi++]=indResDau;
720             if(nFoundKpi>4) return -1;
721           }
722           if(TMath::Abs(pdgresdau)==211){
723             nPions++;
724               sumPxDau+=resdau->Px();
725               sumPyDau+=resdau->Py();
726               sumPzDau+=resdau->Pz();
727               arrayDauLab[nFoundKpi++]=indResDau;
728               if(nFoundKpi>4) return -1;
729           }
730         }
731       }else{
732         return -1;
733       }
734     }
735     if(nPions!=3) return -1;
736     if(nKaons!=1) return -1;
737     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
738     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
739     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
740     return 2;
741   }
742   return -1;
743 }
744
745 //____________________________________________________________________________
746 Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
747  // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
748
749   Int_t pdgD=mcPart->GetPdgCode();
750   if(TMath::Abs(pdgD)!=421) return -1;
751
752   Int_t nDau=mcPart->GetNDaughters();
753
754   if(nDau==2){
755     Int_t daughter0 = mcPart->GetDaughter(0);
756     Int_t daughter1 = mcPart->GetDaughter(1);
757     AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
758     AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
759     if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
760     arrayDauLab[0]=daughter0;
761     arrayDauLab[1]=daughter1;
762     Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
763     Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
764     if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
765        !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
766       return -1;
767     }
768     if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
769     if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
770     if((pdgCode0*pdgCode1)>0) return -1;
771     Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
772     Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
773     Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
774     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
775     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
776     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
777     return 1;
778   }
779
780   if(nDau==3 || nDau==4){
781     Int_t nKaons=0;
782     Int_t nPions=0;
783     Double_t sumPxDau=0.;
784     Double_t sumPyDau=0.;
785     Double_t sumPzDau=0.;
786     Int_t nFoundKpi=0;
787     Int_t labelFirstDau = mcPart->GetDaughter(0);
788     for(Int_t iDau=0; iDau<nDau; iDau++){
789       Int_t indDau = labelFirstDau+iDau;
790       if(indDau<0) return -1;
791       AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
792       if(!dau) return -1;
793       Int_t pdgdau=dau->GetPdgCode();
794       if(TMath::Abs(pdgdau)==321){
795         if(pdgD>0 && pdgdau>0) return -1;
796         if(pdgD<0 && pdgdau<0) return -1;
797         nKaons++;
798         sumPxDau+=dau->Px();
799         sumPyDau+=dau->Py();
800         sumPzDau+=dau->Pz();
801         arrayDauLab[nFoundKpi++]=indDau;
802         if(nFoundKpi>4) return -1;
803       }else if(TMath::Abs(pdgdau)==211){
804         nPions++;
805         sumPxDau+=dau->Px();
806         sumPyDau+=dau->Py();
807         sumPzDau+=dau->Pz();
808         arrayDauLab[nFoundKpi++]=indDau;
809         if(nFoundKpi>4) return -1;
810       }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
811         Int_t nResDau=dau->GetNDaughters();
812         if(nResDau!=2) return -1;
813         Int_t indFirstResDau=dau->GetDaughter(0);
814         for(Int_t resDau=0; resDau<2; resDau++){
815           Int_t indResDau=indFirstResDau+resDau;
816           if(indResDau<0) return -1;
817           AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
818           if(!resdau) return -1;
819           Int_t pdgresdau=resdau->GetPdgCode();
820           if(TMath::Abs(pdgresdau)==321){
821             if(pdgD>0 && pdgresdau>0) return -1;
822             if(pdgD<0 && pdgresdau<0) return -1;
823             nKaons++;
824             sumPxDau+=resdau->Px();
825             sumPyDau+=resdau->Py();
826             sumPzDau+=resdau->Pz();
827             arrayDauLab[nFoundKpi++]=indResDau;
828             if(nFoundKpi>4) return -1;
829           }
830           if(TMath::Abs(pdgresdau)==211){
831             nPions++;
832             sumPxDau+=resdau->Px();
833             sumPyDau+=resdau->Py();
834             sumPzDau+=resdau->Pz();
835             arrayDauLab[nFoundKpi++]=indResDau;
836             if(nFoundKpi>4) return -1;
837           }
838         }
839       }else{
840         return -1;
841       }
842     }
843     if(nPions!=3) return -1;
844     if(nKaons!=1) return -1;
845     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
846     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
847     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
848     return 2;
849   }
850   return -1;
851 }
852 //____________________________________________________________________________
853 Int_t AliVertexingHFUtils::CheckDplusDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
854   // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
855
856   if(label<0) return -1;
857   TParticle* mcPart = stack->Particle(label);
858   Int_t pdgD=mcPart->GetPdgCode();
859   if(TMath::Abs(pdgD)!=411) return -1;
860
861   Int_t nDau=mcPart->GetNDaughters();
862   Int_t labelFirstDau = mcPart->GetDaughter(0);
863   Int_t nKaons=0;
864   Int_t nPions=0;
865   Double_t sumPxDau=0.;
866   Double_t sumPyDau=0.;
867   Double_t sumPzDau=0.;
868   Int_t nFoundKpi=0;
869   
870   if(nDau==3 || nDau==2){
871     for(Int_t iDau=0; iDau<nDau; iDau++){
872       Int_t indDau = labelFirstDau+iDau;
873       if(indDau<0) return -1;
874       TParticle* dau=stack->Particle(indDau);
875       if(!dau) return -1;
876       Int_t pdgdau=dau->GetPdgCode();
877       if(TMath::Abs(pdgdau)==321){
878         if(pdgD*pdgdau>0) return -1;
879         nKaons++;
880         sumPxDau+=dau->Px();
881         sumPyDau+=dau->Py();
882         sumPzDau+=dau->Pz();
883         arrayDauLab[nFoundKpi++]=indDau;
884         if(nFoundKpi>3) return -1;
885       }else if(TMath::Abs(pdgdau)==211){
886         if(pdgD*pdgdau<0) return -1;
887         nPions++;
888         sumPxDau+=dau->Px();
889         sumPyDau+=dau->Py();
890         sumPzDau+=dau->Pz();
891         arrayDauLab[nFoundKpi++]=indDau;
892         if(nFoundKpi>3) return -1;
893       }else if(TMath::Abs(pdgdau)==313){
894         Int_t nResDau=dau->GetNDaughters();
895         if(nResDau!=2) return -1;
896         Int_t indFirstResDau=dau->GetDaughter(0);
897         for(Int_t resDau=0; resDau<2; resDau++){
898           Int_t indResDau=indFirstResDau+resDau;
899           if(indResDau<0) return -1;
900           TParticle* resdau=stack->Particle(indResDau);
901           if(!resdau) return -1;
902           Int_t pdgresdau=resdau->GetPdgCode();
903           if(TMath::Abs(pdgresdau)==321){
904             if(pdgD*pdgresdau>0) return -1;
905             sumPxDau+=resdau->Px();
906             sumPyDau+=resdau->Py();
907             sumPzDau+=resdau->Pz();
908             nKaons++;
909             arrayDauLab[nFoundKpi++]=indResDau;
910             if(nFoundKpi>3) return -1;
911           }
912           if(TMath::Abs(pdgresdau)==211){
913             if(pdgD*pdgresdau<0) return -1;
914             sumPxDau+=resdau->Px();
915             sumPyDau+=resdau->Py();
916             sumPzDau+=resdau->Pz();
917             nPions++;
918             arrayDauLab[nFoundKpi++]=indResDau;
919             if(nFoundKpi>3) return -1;
920           }
921           }
922       }else{
923         return -1;
924       }
925     }
926     if(nPions!=2) return -1;
927     if(nKaons!=1) return -1;
928     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
929     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
930     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
931     if(nDau==3) return 1;
932     else if(nDau==2) return 2;
933   }
934   
935   return -1;
936   
937 }
938
939 //____________________________________________________________________________
940 Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
941   // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
942
943   Int_t pdgD=mcPart->GetPdgCode();
944   if(TMath::Abs(pdgD)!=411) return -1;
945
946   Int_t nDau=mcPart->GetNDaughters();
947   Int_t labelFirstDau = mcPart->GetDaughter(0);
948   Int_t nKaons=0;
949   Int_t nPions=0;
950   Double_t sumPxDau=0.;
951   Double_t sumPyDau=0.;
952   Double_t sumPzDau=0.;
953   Int_t nFoundKpi=0;
954
955   if(nDau==3 || nDau==2){
956     for(Int_t iDau=0; iDau<nDau; iDau++){
957       Int_t indDau = labelFirstDau+iDau;
958       if(indDau<0) return -1;
959       AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
960       if(!dau) return -1;
961       Int_t pdgdau=dau->GetPdgCode();
962       if(TMath::Abs(pdgdau)==321){
963         if(pdgD*pdgdau>0) return -1;
964         nKaons++;
965         sumPxDau+=dau->Px();
966         sumPyDau+=dau->Py();
967         sumPzDau+=dau->Pz();
968         arrayDauLab[nFoundKpi++]=indDau;
969         if(nFoundKpi>3) return -1;
970       }else if(TMath::Abs(pdgdau)==211){
971         if(pdgD*pdgdau<0) return -1;
972         nPions++;
973         sumPxDau+=dau->Px();
974         sumPyDau+=dau->Py();
975         sumPzDau+=dau->Pz();
976         arrayDauLab[nFoundKpi++]=indDau;
977         if(nFoundKpi>3) return -1;
978       }else if(TMath::Abs(pdgdau)==313){
979         Int_t nResDau=dau->GetNDaughters();
980         if(nResDau!=2) return -1;
981         Int_t indFirstResDau=dau->GetDaughter(0);
982         for(Int_t resDau=0; resDau<2; resDau++){
983           Int_t indResDau=indFirstResDau+resDau;
984           if(indResDau<0) return -1;
985           AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
986           if(!resdau) return -1;
987           Int_t pdgresdau=resdau->GetPdgCode();
988           if(TMath::Abs(pdgresdau)==321){
989             if(pdgD*pdgresdau>0) return -1;
990             sumPxDau+=resdau->Px();
991             sumPyDau+=resdau->Py();
992             sumPzDau+=resdau->Pz();
993             nKaons++;
994             arrayDauLab[nFoundKpi++]=indResDau;
995             if(nFoundKpi>3) return -1;
996           }
997           if(TMath::Abs(pdgresdau)==211){
998             if(pdgD*pdgresdau<0) return -1;
999             sumPxDau+=resdau->Px();
1000             sumPyDau+=resdau->Py();
1001             sumPzDau+=resdau->Pz();
1002             nPions++;
1003             arrayDauLab[nFoundKpi++]=indResDau;
1004             if(nFoundKpi>3) return -1;
1005           }
1006         }
1007        }else{
1008         return -1;
1009       }
1010     }
1011     if(nPions!=2) return -1;
1012     if(nKaons!=1) return -1;
1013     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1014     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1015     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1016     if(nDau==3) return 1;
1017     else if(nDau==2) return 2;
1018   }
1019
1020   return -1;
1021
1022 }
1023 //____________________________________________________________________________
1024 Int_t AliVertexingHFUtils::CheckDsDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1025   // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
1026
1027   if(label<0) return -1;
1028   TParticle* mcPart = stack->Particle(label);
1029   Int_t pdgD=mcPart->GetPdgCode();
1030   if(TMath::Abs(pdgD)!=431) return -1;
1031
1032   Int_t nDau=mcPart->GetNDaughters();
1033   Int_t labelFirstDau = mcPart->GetDaughter(0);
1034   Int_t nKaons=0;
1035   Int_t nPions=0;
1036   Double_t sumPxDau=0.;
1037   Double_t sumPyDau=0.;
1038   Double_t sumPzDau=0.;
1039   Int_t nFoundKpi=0;
1040   Bool_t isPhi=kFALSE;
1041   Bool_t isk0st=kFALSE;
1042   
1043   if(nDau==3 || nDau==2){
1044     for(Int_t iDau=0; iDau<nDau; iDau++){
1045       Int_t indDau = labelFirstDau+iDau;
1046       if(indDau<0) return -1;
1047       TParticle* dau=stack->Particle(indDau);
1048       if(!dau) return -1;
1049       Int_t pdgdau=dau->GetPdgCode();
1050       if(TMath::Abs(pdgdau)==321){
1051         nKaons++;
1052         sumPxDau+=dau->Px();
1053         sumPyDau+=dau->Py();
1054         sumPzDau+=dau->Pz();
1055         arrayDauLab[nFoundKpi++]=indDau;
1056         if(nFoundKpi>3) return -1;
1057       }else if(TMath::Abs(pdgdau)==211){
1058         nPions++;
1059         sumPxDau+=dau->Px();
1060         sumPyDau+=dau->Py();
1061         sumPzDau+=dau->Pz();
1062         arrayDauLab[nFoundKpi++]=indDau;
1063         if(nFoundKpi>3) return -1;
1064       }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1065         if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1066         if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1067         Int_t nResDau=dau->GetNDaughters();
1068         if(nResDau!=2) return -1;
1069         Int_t indFirstResDau=dau->GetDaughter(0);
1070         for(Int_t resDau=0; resDau<2; resDau++){
1071           Int_t indResDau=indFirstResDau+resDau;
1072           if(indResDau<0) return -1;
1073           TParticle* resdau=stack->Particle(indResDau);
1074           if(!resdau) return -1;
1075           Int_t pdgresdau=resdau->GetPdgCode();
1076           if(TMath::Abs(pdgresdau)==321){
1077             sumPxDau+=resdau->Px();
1078             sumPyDau+=resdau->Py();
1079             sumPzDau+=resdau->Pz();
1080             nKaons++;
1081             arrayDauLab[nFoundKpi++]=indResDau;
1082             if(nFoundKpi>3) return -1;
1083           }
1084           if(TMath::Abs(pdgresdau)==211){
1085             sumPxDau+=resdau->Px();
1086             sumPyDau+=resdau->Py();
1087             sumPzDau+=resdau->Pz();
1088             nPions++;
1089             arrayDauLab[nFoundKpi++]=indResDau;
1090             if(nFoundKpi>3) return -1;
1091           }
1092           }
1093       }else{
1094         return -1;
1095       }
1096     }
1097     if(nPions!=1) return -1;
1098     if(nKaons!=2) return -1;
1099     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1100     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1101     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1102     if(nDau==3) return 3;
1103     else if(nDau==2){
1104       if(isk0st) return 2;
1105       if(isPhi) return 1;
1106     }
1107   }
1108   
1109   return -1;
1110   
1111 }
1112 //____________________________________________________________________________
1113 Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1114   // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
1115
1116   Int_t pdgD=mcPart->GetPdgCode();
1117   if(TMath::Abs(pdgD)!=431) return -1;
1118
1119   Int_t nDau=mcPart->GetNDaughters();
1120   Int_t labelFirstDau = mcPart->GetDaughter(0);
1121   Int_t nKaons=0;
1122   Int_t nPions=0;
1123   Double_t sumPxDau=0.;
1124   Double_t sumPyDau=0.;
1125   Double_t sumPzDau=0.;
1126   Int_t nFoundKpi=0;
1127   Bool_t isPhi=kFALSE;
1128   Bool_t isk0st=kFALSE;
1129   
1130   if(nDau==3 || nDau==2){
1131     for(Int_t iDau=0; iDau<nDau; iDau++){
1132       Int_t indDau = labelFirstDau+iDau;
1133       if(indDau<0) return -1;
1134       AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1135       if(!dau) return -1;
1136       Int_t pdgdau=dau->GetPdgCode();
1137       if(TMath::Abs(pdgdau)==321){
1138         nKaons++;
1139         sumPxDau+=dau->Px();
1140         sumPyDau+=dau->Py();
1141         sumPzDau+=dau->Pz();
1142         arrayDauLab[nFoundKpi++]=indDau;
1143         if(nFoundKpi>3) return -1;
1144       }else if(TMath::Abs(pdgdau)==211){
1145         nPions++;
1146         sumPxDau+=dau->Px();
1147         sumPyDau+=dau->Py();
1148         sumPzDau+=dau->Pz();
1149         arrayDauLab[nFoundKpi++]=indDau;
1150         if(nFoundKpi>3) return -1;
1151       }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1152         if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1153         if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1154         Int_t nResDau=dau->GetNDaughters();
1155         if(nResDau!=2) return -1;
1156         Int_t indFirstResDau=dau->GetDaughter(0);
1157         for(Int_t resDau=0; resDau<2; resDau++){
1158           Int_t indResDau=indFirstResDau+resDau;
1159           if(indResDau<0) return -1;
1160           AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1161           if(!resdau) return -1;
1162           Int_t pdgresdau=resdau->GetPdgCode();
1163           if(TMath::Abs(pdgresdau)==321){
1164             sumPxDau+=resdau->Px();
1165             sumPyDau+=resdau->Py();
1166             sumPzDau+=resdau->Pz();
1167             nKaons++;
1168             arrayDauLab[nFoundKpi++]=indResDau;
1169             if(nFoundKpi>3) return -1;
1170           }
1171           if(TMath::Abs(pdgresdau)==211){
1172             sumPxDau+=resdau->Px();
1173             sumPyDau+=resdau->Py();
1174             sumPzDau+=resdau->Pz();
1175             nPions++;
1176             arrayDauLab[nFoundKpi++]=indResDau;
1177             if(nFoundKpi>3) return -1;
1178           }
1179           }
1180       }else{
1181         return -1;
1182       }
1183     }
1184     if(nPions!=1) return -1;
1185     if(nKaons!=2) return -1;
1186     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1187     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1188     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1189     if(nDau==3) return 3;
1190     else if(nDau==2){
1191       if(isk0st) return 2;
1192       if(isPhi) return 1;
1193     }
1194   }
1195   
1196   return -1;
1197   
1198 }
1199 //____________________________________________________________________________
1200 Int_t AliVertexingHFUtils::CheckDstarDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1201   // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases
1202
1203   if(label<0) return -1;
1204   TParticle* mcPart = stack->Particle(label);
1205   Int_t pdgD=mcPart->GetPdgCode();
1206   if(TMath::Abs(pdgD)!=413) return -1;
1207
1208   Int_t nDau=mcPart->GetNDaughters();
1209   if(nDau!=2) return -1;
1210
1211   Int_t labelFirstDau = mcPart->GetDaughter(0);
1212   Int_t nKaons=0;
1213   Int_t nPions=0;
1214   Double_t sumPxDau=0.;
1215   Double_t sumPyDau=0.;
1216   Double_t sumPzDau=0.;
1217   Int_t nFoundKpi=0;
1218
1219   for(Int_t iDau=0; iDau<nDau; iDau++){
1220     Int_t indDau = labelFirstDau+iDau;
1221     if(indDau<0) return -1;
1222     TParticle* dau=stack->Particle(indDau);
1223     if(!dau) return -1;
1224     Int_t pdgdau=dau->GetPdgCode();
1225     if(TMath::Abs(pdgdau)==421){
1226       Int_t nResDau=dau->GetNDaughters();
1227       if(nResDau!=2) return -1;
1228       Int_t indFirstResDau=dau->GetDaughter(0);
1229       for(Int_t resDau=0; resDau<2; resDau++){
1230         Int_t indResDau=indFirstResDau+resDau;
1231         if(indResDau<0) return -1;
1232         TParticle* resdau=stack->Particle(indResDau);
1233         if(!resdau) return -1;
1234         Int_t pdgresdau=resdau->GetPdgCode();
1235         if(TMath::Abs(pdgresdau)==321){
1236           if(pdgD*pdgresdau>0) return -1;
1237           sumPxDau+=resdau->Px();
1238           sumPyDau+=resdau->Py();
1239           sumPzDau+=resdau->Pz();
1240           nKaons++;
1241           arrayDauLab[nFoundKpi++]=indResDau;
1242           if(nFoundKpi>3) return -1;
1243         }
1244         if(TMath::Abs(pdgresdau)==211){
1245           if(pdgD*pdgresdau<0) return -1;
1246           sumPxDau+=resdau->Px();
1247           sumPyDau+=resdau->Py();
1248           sumPzDau+=resdau->Pz();
1249           nPions++;
1250           arrayDauLab[nFoundKpi++]=indResDau;
1251           if(nFoundKpi>3) return -1;
1252         }
1253       }
1254     }else if(TMath::Abs(pdgdau)==211){
1255       if(pdgD*pdgdau<0) return -1;
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     }
1263   }
1264
1265   if(nPions!=2) return -1;
1266   if(nKaons!=1) return -1;
1267   if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1268   if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1269   if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1270   return 1;
1271   
1272 }
1273
1274 //____________________________________________________________________________
1275 Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1276   // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases
1277
1278   Int_t pdgD=mcPart->GetPdgCode();
1279   if(TMath::Abs(pdgD)!=413) return -1;
1280
1281   Int_t nDau=mcPart->GetNDaughters();
1282   if(nDau!=2) return -1;
1283
1284   Int_t labelFirstDau = mcPart->GetDaughter(0);
1285   Int_t nKaons=0;
1286   Int_t nPions=0;
1287   Double_t sumPxDau=0.;
1288   Double_t sumPyDau=0.;
1289   Double_t sumPzDau=0.;
1290   Int_t nFoundKpi=0;
1291
1292   for(Int_t iDau=0; iDau<nDau; iDau++){
1293     Int_t indDau = labelFirstDau+iDau;
1294     if(indDau<0) return -1;
1295     AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1296     if(!dau) return -1;
1297     Int_t pdgdau=dau->GetPdgCode();
1298     if(TMath::Abs(pdgdau)==421){
1299       Int_t nResDau=dau->GetNDaughters();
1300       if(nResDau!=2) return -1;
1301       Int_t indFirstResDau=dau->GetDaughter(0);
1302       for(Int_t resDau=0; resDau<2; resDau++){
1303         Int_t indResDau=indFirstResDau+resDau;
1304         if(indResDau<0) return -1;
1305         AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau)); 
1306         if(!resdau) return -1;
1307         Int_t pdgresdau=resdau->GetPdgCode();
1308         if(TMath::Abs(pdgresdau)==321){
1309           if(pdgD*pdgresdau>0) return -1;
1310           sumPxDau+=resdau->Px();
1311           sumPyDau+=resdau->Py();
1312           sumPzDau+=resdau->Pz();
1313           nKaons++;
1314           arrayDauLab[nFoundKpi++]=indResDau;
1315           if(nFoundKpi>3) return -1;
1316         }
1317         if(TMath::Abs(pdgresdau)==211){
1318           if(pdgD*pdgresdau<0) return -1;
1319           sumPxDau+=resdau->Px();
1320           sumPyDau+=resdau->Py();
1321           sumPzDau+=resdau->Pz();
1322           nPions++;
1323           arrayDauLab[nFoundKpi++]=indResDau;
1324           if(nFoundKpi>3) return -1;
1325         }
1326       }
1327     }else if(TMath::Abs(pdgdau)==211){
1328       if(pdgD*pdgdau<0) return -1;
1329       nPions++;
1330       sumPxDau+=dau->Px();
1331       sumPyDau+=dau->Py();
1332       sumPzDau+=dau->Pz();
1333       arrayDauLab[nFoundKpi++]=indDau;
1334       if(nFoundKpi>3) return -1;
1335     }
1336   }
1337
1338   if(nPions!=2) return -1;
1339   if(nKaons!=1) return -1;
1340   if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1341   if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1342   if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1343   return 1;
1344   
1345 }
1346