New methods for tagging decay channels
[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 "AliAODEvent.h"
24 #include "AliAODMCHeader.h"
25 #include "AliGenEventHeader.h"
26 #include "AliAODMCParticle.h"
27 #include "AliAODRecoDecayHF.h"
28 #include "AliVertexingHFUtils.h"
29
30 /* $Id$ */
31
32 ///////////////////////////////////////////////////////////////////
33 //                                                               //
34 // Class with functions useful for different D2H analyses        //
35 // - event plane resolution                                      //
36 // - <pt> calculation with side band subtraction                 //
37 // - tracklet multiplicity calculation                            //
38 // Origin: F.Prino, Torino, prino@to.infn.it                     //
39 //                                                               //
40 ///////////////////////////////////////////////////////////////////
41
42 ClassImp(AliVertexingHFUtils)
43
44 //______________________________________________________________________
45 AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
46   fK(1),
47   fSubRes(1.),
48   fMinEtaForTracklets(-1.),
49   fMaxEtaForTracklets(1.)
50 {
51   // Default contructor
52 }
53
54
55 //______________________________________________________________________
56 AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
57   TObject(),
58   fK(k),
59   fSubRes(1.),
60   fMinEtaForTracklets(-1.),
61   fMaxEtaForTracklets(1.)
62 {
63   // Standard constructor
64 }
65
66
67 //______________________________________________________________________
68 void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t  errsignal, Double_t  background, Double_t  errbackground, Double_t &significance,Double_t &errsignificance){
69   // calculate significance from S, B and errors
70
71
72   Double_t errSigSq=errsignal*errsignal;
73   Double_t errBkgSq=errbackground*errbackground;
74   Double_t sigPlusBkg=signal+background;
75   if(sigPlusBkg>0. && signal>0.){
76     significance =  signal/TMath::Sqrt(signal+background);
77     errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal);
78   }else{
79     significance=0.;
80     errsignificance=0.;
81   }
82   return;
83
84 }
85 //______________________________________________________________________
86 Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
87   // compute chi from polynomial approximation
88   Double_t c[5];
89   if(k==1){ 
90     c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
91   }
92   else if(k==2){
93     c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
94   }
95   else return -1;
96   return c[0]*x+c[1]*x*x+c[2]*x*x*x+c[3]*x*x*x*x+c[4]*x*x*x*x*x;
97 }
98
99 //______________________________________________________________________
100 Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
101   return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
102 }
103
104
105 //______________________________________________________________________
106 Double_t AliVertexingHFUtils::FindChi(Double_t res,  Int_t k){
107   // compute chi variable (=v2*sqrt(N)) from external values
108
109   Double_t x1=0;
110   Double_t x2=15;
111   Double_t y1,y2;
112   if(k==1){
113     y1=ResolK1(x1)-res;
114     y2=ResolK1(x2)-res;
115   }
116   else if(k==2){
117     y1=Pol(x1,2)-res;
118     y2=Pol(x2,2)-res;
119   }
120   else return -1;
121
122   if(y1*y2>0) return -1;
123   if(y1==0) return y1;
124   if(y2==0) return y2;
125   Double_t xmed,ymed;
126   Int_t jiter=0;
127   while((x2-x1)>0.0001){
128     xmed=0.5*(x1+x2);
129     if(k==1){
130       y1=ResolK1(x1)-res;
131       y2=ResolK1(x2)-res;
132       ymed=ResolK1(xmed)-res;
133     }
134     else if(k==2){
135       y1=Pol(x1,2)-res;
136       y2=Pol(x2,2)-res;
137       ymed=Pol(xmed,2)-res;
138     }
139     else return -1;
140     if((y1*ymed)<0) x2=xmed;
141     if((y2*ymed)<0)x1=xmed;
142     if(ymed==0) return xmed;
143     jiter++;
144   }
145   return 0.5*(x1+x2);
146 }
147
148 //______________________________________________________________________
149 Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
150   // computes event plane resolution starting from sub event resolution
151   Double_t chisub=FindChi(resSub,k);
152   Double_t chifull=chisub*TMath::Sqrt(2);
153   if(k==1) return ResolK1(chifull);
154   else if(k==2) return Pol(chifull,2);
155   else{
156     printf("k should be <=2\n");
157     return 1.;
158   }
159 }
160
161 //______________________________________________________________________
162 Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
163   // computes event plane resolution starting from sub event correlation histogram
164   if(!hSubEvCorr) return 1.;
165   Double_t resSub=GetSubEvResol(hSubEvCorr);
166   return GetFullEvResol(resSub,k);
167 }
168 //______________________________________________________________________
169 Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
170   // computes low limit event plane resolution starting from sub event correlation histogram
171   if(!hSubEvCorr) return 1.;
172   Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
173   printf("%f\n",resSub);
174   return GetFullEvResol(resSub,k);  
175 }
176 //______________________________________________________________________
177 Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
178   // computes high limit event plane resolution starting from sub event correlation histogram
179   if(!hSubEvCorr) return 1.;
180   Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
181   printf("%f\n",resSub);
182   return GetFullEvResol(resSub,k);  
183 }
184 //______________________________________________________________________
185 Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
186   // counts tracklets in given eta range
187   AliAODTracklets* tracklets=ev->GetTracklets();
188   Int_t nTr=tracklets->GetNumberOfTracklets();
189   Int_t count=0;
190   for(Int_t iTr=0; iTr<nTr; iTr++){
191     Double_t theta=tracklets->GetTheta(iTr);
192     Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
193     if(eta>mineta && eta<maxeta) count++;
194   }
195   return count;
196 }
197 //______________________________________________________________________
198 Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
199   // counts generated particles in fgiven eta range
200
201   Int_t nChargedMC=0;
202   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
203     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
204     Int_t charge = part->Charge();
205     Double_t eta = part->Eta();
206     if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
207   } 
208   return nChargedMC;
209 }
210 //______________________________________________________________________
211 Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
212   // counts generated primary particles in given eta range
213
214   Int_t nChargedMC=0;
215   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
216     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
217     Int_t charge = part->Charge();
218     Double_t eta = part->Eta();
219     if(charge!=0 && eta>mineta && eta<maxeta){
220       if(part->IsPrimary())nChargedMC++;
221     } 
222   }  
223   return nChargedMC;
224 }
225 //______________________________________________________________________
226 Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
227   // counts generated primary particles in given eta range
228
229   Int_t nChargedMC=0;
230   for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
231     AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
232     Int_t charge = part->Charge();
233     Double_t eta = part->Eta();
234     if(charge!=0 && eta>mineta && eta<maxeta){
235       if(part->IsPhysicalPrimary())nChargedMC++;
236     } 
237   }
238   return nChargedMC;
239 }
240 //______________________________________________________________________
241 void AliVertexingHFUtils::AveragePt(Float_t& averagePt, Float_t& errorPt,Float_t ptmin,Float_t ptmax, TH2F* hMassD, Float_t massFromFit, Float_t sigmaFromFit, TF1* funcB2, Float_t sigmaRangeForSig,Float_t sigmaRangeForBkg, Float_t minMass, Float_t maxMass, Int_t rebin){
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;
266   Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2);
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);
272   Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1);
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
316 }
317 //____________________________________________________________________________
318 Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
319   // true impact parameter calculation for Dzero
320
321   if(!partD || !arrayMC || !mcHeader) return 99999.;
322   Int_t code=partD->GetPdgCode();
323   if(TMath::Abs(code)!=421) return 99999.;
324
325   Double_t vtxTrue[3];
326   mcHeader->GetVertex(vtxTrue);
327   Double_t origD[3];
328   partD->XvYvZv(origD);
329   Short_t charge=partD->Charge();
330   Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
331   for(Int_t iDau=0; iDau<2; iDau++){
332     pXdauTrue[iDau]=0.;
333     pYdauTrue[iDau]=0.;
334     pZdauTrue[iDau]=0.;
335   }
336
337   Int_t nDau=partD->GetNDaughters();
338   Int_t labelFirstDau = partD->GetDaughter(0); 
339   if(nDau==2){
340     for(Int_t iDau=0; iDau<2; iDau++){
341       Int_t ind = labelFirstDau+iDau;
342       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
343       if(!part){
344         printf("Daughter particle not found in MC array");
345         return 99999.;
346       } 
347       pXdauTrue[iDau]=part->Px();
348       pYdauTrue[iDau]=part->Py();
349       pZdauTrue[iDau]=part->Pz();
350     }
351   }else{
352     return 99999.;
353   }
354
355   Double_t d0dummy[2]={0.,0.};
356   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
357   return aodDvsMC.ImpParXY();
358
359 }
360 //____________________________________________________________________________
361 Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
362   // true impact parameter calculation for Dplus
363
364   if(!partD || !arrayMC || !mcHeader) return 99999.;
365   Int_t code=partD->GetPdgCode();
366   if(TMath::Abs(code)!=411) return 99999.;
367
368   Double_t vtxTrue[3];
369   mcHeader->GetVertex(vtxTrue);
370   Double_t origD[3];
371   partD->XvYvZv(origD);
372   Short_t charge=partD->Charge();
373   Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
374   for(Int_t iDau=0; iDau<3; iDau++){
375     pXdauTrue[iDau]=0.;
376     pYdauTrue[iDau]=0.;
377     pZdauTrue[iDau]=0.;
378   }
379
380   Int_t nDau=partD->GetNDaughters();
381   Int_t labelFirstDau = partD->GetDaughter(0); 
382   if(nDau==3){
383     for(Int_t iDau=0; iDau<3; iDau++){
384       Int_t ind = labelFirstDau+iDau;
385       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
386       if(!part){
387         printf("Daughter particle not found in MC array");
388         return 99999.;
389       } 
390       pXdauTrue[iDau]=part->Px();
391       pYdauTrue[iDau]=part->Py();
392       pZdauTrue[iDau]=part->Pz();
393     }
394   }else if(nDau==2){
395     Int_t theDau=0;
396     for(Int_t iDau=0; iDau<2; iDau++){
397       Int_t ind = labelFirstDau+iDau;
398       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
399       if(!part){
400         printf("Daughter particle not found in MC array");
401         return 99999.;
402       } 
403       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
404       if(pdgCode==211 || pdgCode==321){
405         pXdauTrue[theDau]=part->Px();
406         pYdauTrue[theDau]=part->Py();
407         pZdauTrue[theDau]=part->Pz();
408         ++theDau;
409       }else{
410         Int_t nDauRes=part->GetNDaughters();
411         if(nDauRes==2){
412           Int_t labelFirstDauRes = part->GetDaughter(0);        
413           for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
414             Int_t indDR = labelFirstDauRes+iDauRes;
415             AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
416             if(!partDR){
417               printf("Daughter particle not found in MC array");
418               return 99999.;
419             } 
420             
421             Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
422             if(pdgCodeDR==211 || pdgCodeDR==321){
423               pXdauTrue[theDau]=partDR->Px();
424               pYdauTrue[theDau]=partDR->Py();
425               pZdauTrue[theDau]=partDR->Pz();
426               ++theDau;
427             }
428           }
429         }
430       }
431     }
432     if(theDau!=3){
433       printf("Wrong number of decay prongs");
434       return 99999.;
435     }
436   }
437
438   Double_t d0dummy[3]={0.,0.,0.};
439   AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
440   return aodDvsMC.ImpParXY();
441
442 }
443
444
445
446 //____________________________________________________________________________
447 Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
448   //
449   // Correct the number of accepted tracklets based on a TProfile Hist
450   //
451   //
452
453   if(TMath::Abs(vtxZ)>10.0){
454     //    printf("ERROR: Z vertex out of range for correction of multiplicity\n");
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 }
471 //______________________________________________________________________
472 TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
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 }
488 //_____________________________________________________________________
489 void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
490
491   // method to check if a track comes from a given generator
492
493   Int_t lab=TMath::Abs(track->GetLabel());
494   nameGen=GetGenerator(lab,header);
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   
518   return;
519 }
520 //----------------------------------------------------------------------
521 Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
522   // method to check if a track comes from the signal event or from the underlying Hijing event
523   TString nameGen;
524   GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
525   
526   if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
527   
528   return kTRUE;
529 }
530 //____________________________________________________________________________
531 Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
532   // method to check if a D meson candidate comes from the signal event or from the underlying Hijing event
533
534   Int_t nprongs=cand->GetNProngs();
535   for(Int_t i=0;i<nprongs;i++){
536     AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i);
537     if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
538   }
539   return kFALSE;
540 }
541 //____________________________________________________________________________
542 Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
543   // checking whether the mother of the particles come from a charm or a bottom quark
544
545   Int_t pdgGranma = 0;
546   Int_t mother = 0;
547   mother = mcPart->GetMother();
548   Int_t istep = 0;
549   Int_t abspdgGranma =0;
550   Bool_t isFromB=kFALSE;
551   Bool_t isQuarkFound=kFALSE;
552   while (mother >0 ){
553     istep++;
554     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
555     if (mcGranma){
556       pdgGranma = mcGranma->GetPdgCode();
557       abspdgGranma = TMath::Abs(pdgGranma);
558       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
559         isFromB=kTRUE;
560       }
561       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
562       mother = mcGranma->GetMother();
563     }else{
564       printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
565       break;
566     }
567   }
568   if(searchUpToQuark && !isQuarkFound) return 0;
569   if(isFromB) return 5;
570   else return 4;
571
572 }
573 //____________________________________________________________________________
574 Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
575   // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
576
577
578   Int_t pdgD=mcPart->GetPdgCode();
579   if(TMath::Abs(pdgD)!=421) return -1;
580
581   Int_t nDau=mcPart->GetNDaughters();
582
583   if(nDau==2){
584     Int_t daughter0 = mcPart->GetDaughter(0);
585     Int_t daughter1 = mcPart->GetDaughter(1);
586     AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
587     AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
588     if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
589     arrayDauLab[0]=daughter0;
590     arrayDauLab[1]=daughter1;
591     Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
592     Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
593     if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
594        !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
595       return -1;
596     }
597     if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
598     if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
599     if((pdgCode0*pdgCode1)>0) return -1;
600     Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
601     Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
602     Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
603     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
604     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
605     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
606     return 1;
607   }
608
609   if(nDau==3 || nDau==4){
610     Int_t nKaons=0;
611     Int_t nPions=0;
612     Double_t sumPxDau=0.;
613     Double_t sumPyDau=0.;
614     Double_t sumPzDau=0.;
615     Int_t nFoundKpi=0;
616     Int_t labelFirstDau = mcPart->GetDaughter(0);
617     for(Int_t iDau=0; iDau<nDau; iDau++){
618       Int_t indDau = labelFirstDau+iDau;
619       if(indDau<0) return -1;
620       AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
621       if(!dau) return -1;
622       Int_t pdgdau=dau->GetPdgCode();
623       if(TMath::Abs(pdgdau)==321){
624         if(pdgD>0 && pdgdau>0) return -1;
625         if(pdgD<0 && pdgdau<0) return -1;
626         nKaons++;
627         sumPxDau+=dau->Px();
628         sumPyDau+=dau->Py();
629         sumPzDau+=dau->Pz();
630         arrayDauLab[nFoundKpi++]=indDau;
631         if(nFoundKpi>4) return -1;
632       }else if(TMath::Abs(pdgdau)==211){
633         nPions++;
634         sumPxDau+=dau->Px();
635         sumPyDau+=dau->Py();
636         sumPzDau+=dau->Pz();
637         arrayDauLab[nFoundKpi++]=indDau;
638         if(nFoundKpi>4) return -1;
639       }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
640         Int_t nResDau=dau->GetNDaughters();
641         if(nResDau!=2) return -1;
642         Int_t indFirstResDau=dau->GetDaughter(0);
643         for(Int_t resDau=0; resDau<2; resDau++){
644           Int_t indResDau=indFirstResDau+resDau;
645           if(indResDau<0) return -1;
646           AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
647           if(!resdau) return -1;
648           Int_t pdgresdau=resdau->GetPdgCode();
649           if(TMath::Abs(pdgresdau)==321){
650             if(pdgD>0 && pdgresdau>0) return -1;
651             if(pdgD<0 && pdgresdau<0) return -1;
652             nKaons++;
653             sumPxDau+=resdau->Px();
654             sumPyDau+=resdau->Py();
655             sumPzDau+=resdau->Pz();
656             arrayDauLab[nFoundKpi++]=indResDau;
657             if(nFoundKpi>4) return -1;
658           }
659           if(TMath::Abs(pdgresdau)==211){
660             nPions++;
661             sumPxDau+=resdau->Px();
662             sumPyDau+=resdau->Py();
663             sumPzDau+=resdau->Pz();
664             arrayDauLab[nFoundKpi++]=indResDau;
665             if(nFoundKpi>4) return -1;
666           }
667         }
668       }else{
669         return -1;
670       }
671     }
672     if(nPions!=3) return -1;
673     if(nKaons!=1) return -1;
674     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
675     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
676     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
677     return 2;
678   }
679   return -1;
680 }
681 //____________________________________________________________________________
682 Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
683   // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
684
685   Int_t pdgD=mcPart->GetPdgCode();
686   if(TMath::Abs(pdgD)!=411) return -1;
687
688   Int_t nDau=mcPart->GetNDaughters();
689   Int_t labelFirstDau = mcPart->GetDaughter(0);
690   Int_t nKaons=0;
691   Int_t nPions=0;
692   Double_t sumPxDau=0.;
693   Double_t sumPyDau=0.;
694   Double_t sumPzDau=0.;
695   Int_t nFoundKpi=0;
696
697   if(nDau==3 || nDau==2){
698     for(Int_t iDau=0; iDau<nDau; iDau++){
699       Int_t indDau = labelFirstDau+iDau;
700       if(indDau<0) return -1;
701       AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
702       if(!dau) return -1;
703       Int_t pdgdau=dau->GetPdgCode();
704       if(TMath::Abs(pdgdau)==321){
705         if(pdgD*pdgdau>0) return -1;
706         nKaons++;
707         sumPxDau+=dau->Px();
708         sumPyDau+=dau->Py();
709         sumPzDau+=dau->Pz();
710         arrayDauLab[nFoundKpi++]=indDau;
711         if(nFoundKpi>3) return -1;
712       }else if(TMath::Abs(pdgdau)==211){
713         if(pdgD*pdgdau<0) return -1;
714         nPions++;
715         sumPxDau+=dau->Px();
716         sumPyDau+=dau->Py();
717         sumPzDau+=dau->Pz();
718         arrayDauLab[nFoundKpi++]=indDau;
719         if(nFoundKpi>3) return -1;
720       }else if(TMath::Abs(pdgdau)==313){
721         Int_t nResDau=dau->GetNDaughters();
722         if(nResDau!=2) return -1;
723         Int_t indFirstResDau=dau->GetDaughter(0);
724         for(Int_t resDau=0; resDau<2; resDau++){
725           Int_t indResDau=indFirstResDau+resDau;
726           if(indResDau<0) return -1;
727           AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
728           if(!resdau) return -1;
729           Int_t pdgresdau=resdau->GetPdgCode();
730           if(TMath::Abs(pdgresdau)==321){
731             if(pdgD*pdgresdau>0) return -1;
732             sumPxDau+=resdau->Px();
733             sumPyDau+=resdau->Py();
734             sumPzDau+=resdau->Pz();
735             nKaons++;
736             arrayDauLab[nFoundKpi++]=indResDau;
737             if(nFoundKpi>3) return -1;
738           }
739           if(TMath::Abs(pdgresdau)==211){
740             if(pdgD*pdgresdau<0) return -1;
741             sumPxDau+=resdau->Px();
742             sumPyDau+=resdau->Py();
743             sumPzDau+=resdau->Pz();
744             nPions++;
745             arrayDauLab[nFoundKpi++]=indResDau;
746             if(nFoundKpi>3) return -1;
747           }
748         }
749        }else{
750         return -1;
751       }
752     }
753     if(nPions!=2) return -1;
754     if(nKaons!=1) return -1;
755     if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
756     if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
757     if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
758     if(nDau==3) return 1;
759     else if(nDau==2) return 2;
760   }
761
762   return -1;
763
764 }