AliPhysics  608b256 (608b256)
AliVertexingHFUtils.cxx
Go to the documentation of this file.
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 "AliMCEvent.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 
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 // //
43 
45 ClassImp(AliVertexingHFUtils);
47 
48 using namespace std;
49 
50 //______________________________________________________________________
52  fK(1),
53  fSubRes(1.),
54  fMinEtaForTracklets(-1.),
55  fMaxEtaForTracklets(1.)
56 {
58 }
59 
60 
61 //______________________________________________________________________
63  TObject(),
64  fK(k),
65  fSubRes(1.),
68 {
70 }
71 
72 
73 //______________________________________________________________________
74 void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance,Double_t &errsignificance){
76 
77 
78  Double_t errSigSq=errsignal*errsignal;
79  Double_t errBkgSq=errbackground*errbackground;
80  Double_t sigPlusBkg=signal+background;
81  if(sigPlusBkg>0. && signal>0.){
82  significance = signal/TMath::Sqrt(signal+background);
83  errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal);
84  }else{
85  significance=0.;
86  errsignificance=0.;
87  }
88  return;
89 
90 }
91 //______________________________________________________________________
94  Double_t c[5];
95  if(k==1){
96  c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
97  }
98  else if(k==2){
99  c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
100  }
101  else return -1;
102  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;
103 }
104 
105 //______________________________________________________________________
107  return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
108 }
109 
110 
111 //______________________________________________________________________
114 
115  Double_t x1=0;
116  Double_t x2=15;
117  Double_t y1,y2;
118  if(k==1){
119  y1=ResolK1(x1)-res;
120  y2=ResolK1(x2)-res;
121  }
122  else if(k==2){
123  y1=Pol(x1,2)-res;
124  y2=Pol(x2,2)-res;
125  }
126  else return -1;
127 
128  if(y1*y2>0) return -1;
129  if(y1==0) return y1;
130  if(y2==0) return y2;
131  Double_t xmed,ymed;
132  Int_t jiter=0;
133  while((x2-x1)>0.0001){
134  xmed=0.5*(x1+x2);
135  if(k==1){
136  y1=ResolK1(x1)-res;
137  y2=ResolK1(x2)-res;
138  ymed=ResolK1(xmed)-res;
139  }
140  else if(k==2){
141  y1=Pol(x1,2)-res;
142  y2=Pol(x2,2)-res;
143  ymed=Pol(xmed,2)-res;
144  }
145  else return -1;
146  if((y1*ymed)<0) x2=xmed;
147  if((y2*ymed)<0)x1=xmed;
148  if(ymed==0) return xmed;
149  jiter++;
150  }
151  return 0.5*(x1+x2);
152 }
153 
154 //______________________________________________________________________
157  Double_t chisub=FindChi(resSub,k);
158  Double_t chifull=chisub*TMath::Sqrt(2);
159  if(k==1) return ResolK1(chifull);
160  else if(k==2) return Pol(chifull,2);
161  else{
162  printf("k should be <=2\n");
163  return 1.;
164  }
165 }
166 
167 //______________________________________________________________________
170  if(!hSubEvCorr) return 1.;
171  Double_t resSub=GetSubEvResol(hSubEvCorr);
172  return GetFullEvResol(resSub,k);
173 }
174 //______________________________________________________________________
177  if(!hSubEvCorr) return 1.;
178  Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
179  printf("%f\n",resSub);
180  return GetFullEvResol(resSub,k);
181 }
182 //______________________________________________________________________
185  if(!hSubEvCorr) return 1.;
186  Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
187  printf("%f\n",resSub);
188  return GetFullEvResol(resSub,k);
189 }
190 //______________________________________________________________________
193  AliAODTracklets* tracklets=ev->GetTracklets();
194  Int_t nTr=tracklets->GetNumberOfTracklets();
195  Int_t count=0;
196  for(Int_t iTr=0; iTr<nTr; iTr++){
197  Double_t theta=tracklets->GetTheta(iTr);
198  Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
199  if(eta>mineta && eta<maxeta) count++;
200  }
201  return count;
202 }
203 //______________________________________________________________________
206 
207  Int_t nChargedMC=0;
208  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
209  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
210  Int_t charge = part->Charge();
211  Double_t eta = part->Eta();
212  if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
213  }
214  return nChargedMC;
215 }
216 //______________________________________________________________________
219 
220  Int_t nChargedMC=0;
221  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
222  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
223  Int_t charge = part->Charge();
224  Double_t eta = part->Eta();
225  if(charge!=0 && eta>mineta && eta<maxeta){
226  if(part->IsPrimary())nChargedMC++;
227  }
228  }
229  return nChargedMC;
230 }
231 //______________________________________________________________________
234 
235  Int_t nChargedMC=0;
236  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
237  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
238  Int_t charge = part->Charge();
239  Double_t eta = part->Eta();
240  if(charge!=0 && eta>mineta && eta<maxeta){
241  if(part->IsPhysicalPrimary())nChargedMC++;
242  }
243  }
244  return nChargedMC;
245 }
246 
247 
248 //______________________________________________________________________
250  //
253  // http://git.cern.ch/pubweb/AliRoot.git/blob/HEAD:/ANALYSIS/AliCentralitySelectionTask.cxx#l1345
254 
255  Double_t multV0AEq=0;
256  for(Int_t iCh = 32; iCh < 64; ++iCh) {
257  Double_t mult = ev->GetVZEROEqMultiplicity(iCh);
258  multV0AEq += mult;
259  }
260  return multV0AEq;
261 }
262 
263 //______________________________________________________________________
268 
269  Double_t multV0CEq=0;
270  for(Int_t iCh = 0; iCh < 32; ++iCh) {
271  Double_t mult = ev->GetVZEROEqMultiplicity(iCh);
272  multV0CEq += mult;
273  }
274  return multV0CEq;
275 }
276 
277 //______________________________________________________________________
278 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){
279 
281 
282  //Make 2D histos in the desired pt range
283  Int_t start=hMassD->FindBin(ptmin);
284  Int_t end=hMassD->FindBin(ptmax)-1;
285  const Int_t nx=end-start;
286  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()));
287  for(Int_t ix=start;ix<end;ix++){
288  for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
289  hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
290  hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
291  }
292  }
293 
294  Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
295  Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
296  Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
297  Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
298  Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
299  Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
300  // printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
301 
302  Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
303  Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2);
304  Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
305  Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
306  Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
307  Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
308  Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
309  Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1);
310  Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
311  Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
312  // printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
313 
314  Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
315  Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
316  Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
317  // printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
318 
319  TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
320  TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
321  TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
322 
323  hMptBkgLo->Rebin(rebin);
324  hMptBkgHi->Rebin(rebin);
325  hMptSigReg->Rebin(rebin);
326 
327  hMptBkgLo->Sumw2();
328  hMptBkgHi->Sumw2();
329  TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
330  hMptBkgLoScal->Scale(bkgSig/bkgLow);
331  TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
332  hMptBkgHiScal->Scale(bkgSig/bkgHi);
333 
334  TH1F* hMptBkgAver=0x0;
335  hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
336  hMptBkgAver->Add(hMptBkgHiScal);
337  hMptBkgAver->Scale(0.5);
338  TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
339  hMptSig->Add(hMptBkgAver,-1.);
340 
341  averagePt = hMptSig->GetMean();
342  errorPt = hMptSig->GetMeanError();
343 
344  delete hMptBkgLo;
345  delete hMptBkgHi;
346  delete hMptSigReg;
347  delete hMptBkgLoScal;
348  delete hMptBkgHiScal;
349  delete hMptBkgAver;
350  delete hMassDpt;
351  delete hMptSig;
352 
353 }
354 //____________________________________________________________________________
357  const Double32_t *mean = aodEv->GetT0TOF();
358  if(mean && mean[0]<9999.) return kTRUE;
359  else return kFALSE;
360 }
361 //____________________________________________________________________________
362 Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
364 
365  if(!partD || !arrayMC || !mcHeader) return 99999.;
366  Int_t code=partD->GetPdgCode();
367  if(TMath::Abs(code)!=421) return 99999.;
368 
369  Double_t vtxTrue[3];
370  mcHeader->GetVertex(vtxTrue);
371  Double_t origD[3];
372  partD->XvYvZv(origD);
373  Short_t charge=partD->Charge();
374  Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
375  for(Int_t iDau=0; iDau<2; iDau++){
376  pXdauTrue[iDau]=0.;
377  pYdauTrue[iDau]=0.;
378  pZdauTrue[iDau]=0.;
379  }
380 
381  Int_t nDau=partD->GetNDaughters();
382  Int_t labelFirstDau = partD->GetDaughterLabel(0);
383  if(nDau==2){
384  for(Int_t iDau=0; iDau<2; iDau++){
385  Int_t ind = labelFirstDau+iDau;
386  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
387  if(!part){
388  printf("Daughter particle not found in MC array");
389  return 99999.;
390  }
391  pXdauTrue[iDau]=part->Px();
392  pYdauTrue[iDau]=part->Py();
393  pZdauTrue[iDau]=part->Pz();
394  }
395  }else{
396  return 99999.;
397  }
398 
399  Double_t d0dummy[2]={0.,0.};
400  AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
401  return aodDvsMC.ImpParXY();
402 
403 }
404 //____________________________________________________________________________
405 Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
407 
408  if(!partD || !arrayMC || !mcHeader) return 99999.;
409  Int_t code=partD->GetPdgCode();
410  if(TMath::Abs(code)!=411) return 99999.;
411 
412  Double_t vtxTrue[3];
413  mcHeader->GetVertex(vtxTrue);
414  Double_t origD[3];
415  partD->XvYvZv(origD);
416  Short_t charge=partD->Charge();
417  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
418  for(Int_t iDau=0; iDau<3; iDau++){
419  pXdauTrue[iDau]=0.;
420  pYdauTrue[iDau]=0.;
421  pZdauTrue[iDau]=0.;
422  }
423 
424  Int_t nDau=partD->GetNDaughters();
425  Int_t labelFirstDau = partD->GetDaughterLabel(0);
426  if(nDau==3){
427  for(Int_t iDau=0; iDau<3; iDau++){
428  Int_t ind = labelFirstDau+iDau;
429  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
430  if(!part){
431  printf("Daughter particle not found in MC array");
432  return 99999.;
433  }
434  pXdauTrue[iDau]=part->Px();
435  pYdauTrue[iDau]=part->Py();
436  pZdauTrue[iDau]=part->Pz();
437  }
438  }else if(nDau==2){
439  Int_t theDau=0;
440  for(Int_t iDau=0; iDau<2; iDau++){
441  Int_t ind = labelFirstDau+iDau;
442  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
443  if(!part){
444  printf("Daughter particle not found in MC array");
445  return 99999.;
446  }
447  Int_t pdgCode=TMath::Abs(part->GetPdgCode());
448  if(pdgCode==211 || pdgCode==321){
449  pXdauTrue[theDau]=part->Px();
450  pYdauTrue[theDau]=part->Py();
451  pZdauTrue[theDau]=part->Pz();
452  ++theDau;
453  }else{
454  Int_t nDauRes=part->GetNDaughters();
455  if(nDauRes==2){
456  Int_t labelFirstDauRes = part->GetDaughterLabel(0);
457  for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
458  Int_t indDR = labelFirstDauRes+iDauRes;
459  AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
460  if(!partDR){
461  printf("Daughter particle not found in MC array");
462  return 99999.;
463  }
464 
465  Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
466  if(pdgCodeDR==211 || pdgCodeDR==321){
467  pXdauTrue[theDau]=partDR->Px();
468  pYdauTrue[theDau]=partDR->Py();
469  pZdauTrue[theDau]=partDR->Pz();
470  ++theDau;
471  }
472  }
473  }
474  }
475  }
476  if(theDau!=3){
477  printf("Wrong number of decay prongs");
478  return 99999.;
479  }
480  }
481 
482  Double_t d0dummy[3]={0.,0.,0.};
483  AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
484  return aodDvsMC.ImpParXY();
485 
486 }
487 
488 
489 
490 //____________________________________________________________________________
491 Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
492  //
493  // Correct the number of accepted tracklets based on a TProfile Hist
494  //
495  //
496 
497  if(TMath::Abs(vtxZ)>10.0){
498  // printf("ERROR: Z vertex out of range for correction of multiplicity\n");
499  return uncorrectedNacc;
500  }
501 
502  if(!estimatorAvg){
503  printf("ERROR: Missing TProfile for correction of multiplicity\n");
504  return uncorrectedNacc;
505  }
506 
507  Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
508 
509  Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
510 
511  Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
512 
513  return correctedNacc;
514 }
515 //______________________________________________________________________
516 TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
518 
519  Int_t nsumpart=0;
520  TList *lh=header->GetCocktailHeaders();
521  Int_t nh=lh->GetEntries();
522  for(Int_t i=0;i<nh;i++){
523  AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
524  TString genname=gh->GetName();
525  Int_t npart=gh->NProduced();
526  if(label>=nsumpart && label<(nsumpart+npart)) return genname;
527  nsumpart+=npart;
528  }
529  TString empty="";
530  return empty;
531 }
532 //_____________________________________________________________________
533 void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
534 
536 
537  Int_t lab=TMath::Abs(track->GetLabel());
538  nameGen=GetGenerator(lab,header);
539 
540  // Int_t countControl=0;
541 
542  while(nameGen.IsWhitespace()){
543  AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
544  if(!mcpart){
545  printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
546  break;
547  }
548  Int_t mother = mcpart->GetMother();
549  if(mother<0){
550  printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
551  break;
552  }
553  lab=mother;
554  nameGen=GetGenerator(mother,header);
555  // countControl++;
556  // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
557  // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
558  // break;
559  // }
560  }
561 
562  return;
563 }
564 //----------------------------------------------------------------------
565 Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
567  TString nameGen;
568  GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
569 
570  if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
571 
572  return kTRUE;
573 }
574 //____________________________________________________________________________
575 Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
577 
578  Int_t nprongs=cand->GetNProngs();
579  for(Int_t i=0;i<nprongs;i++){
580  AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i);
581  if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
582  }
583  return kFALSE;
584 }
585 //____________________________________________________________________________
586 Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
588 
589  AliAODTrack* bach = cand->GetBachelor();
590  if(IsTrackInjected(bach, header, arrayMC)) {
591  AliDebug(2, "Bachelor is injected, the whole candidate is then injected");
592  return kTRUE;
593  }
594  AliAODv0* v0 = cand->Getv0();
595  Int_t nprongs = v0->GetNProngs();
596  for(Int_t i = 0; i < nprongs; i++){
597  AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i);
598  if(IsTrackInjected(daugh,header,arrayMC)) {
599  AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i));
600  return kTRUE;
601  }
602  }
603  return kFALSE;
604 }
605 //____________________________________________________________________________
606 Int_t AliVertexingHFUtils::CheckOrigin(AliMCEvent* mcEvent, AliMCParticle *mcPart, Bool_t searchUpToQuark){
608 
609  Int_t pdgGranma = 0;
610  Int_t mother = 0;
611  mother = mcPart->GetMother();
612  Int_t istep = 0;
613  Int_t abspdgGranma =0;
614  Bool_t isFromB=kFALSE;
615  Bool_t isQuarkFound=kFALSE;
616  while (mother >0 ){
617  istep++;
618  AliMCParticle* mcGranma = (AliMCParticle*)mcEvent->GetTrack(mother);
619  if (mcGranma){
620  TParticle* partGranma = mcGranma->Particle();
621  if(partGranma) pdgGranma = partGranma->GetPdgCode();
622  abspdgGranma = TMath::Abs(pdgGranma);
623  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
624  isFromB=kTRUE;
625  }
626  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
627  mother = mcGranma->GetMother();
628  }else{
629  printf("CheckOrigin: Failed casting the mother particle!");
630  break;
631  }
632  }
633  if(searchUpToQuark && !isQuarkFound) return 0;
634  if(isFromB) return 5;
635  else return 4;
636 
637 }
638 //____________________________________________________________________________
639 Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
641 
642  Int_t pdgGranma = 0;
643  Int_t mother = 0;
644  mother = mcPart->GetMother();
645  Int_t istep = 0;
646  Int_t abspdgGranma =0;
647  Bool_t isFromB=kFALSE;
648  Bool_t isQuarkFound=kFALSE;
649  while (mother >0 ){
650  istep++;
651  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
652  if (mcGranma){
653  pdgGranma = mcGranma->GetPdgCode();
654  abspdgGranma = TMath::Abs(pdgGranma);
655  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
656  isFromB=kTRUE;
657  }
658  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
659  mother = mcGranma->GetMother();
660  }else{
661  printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
662  break;
663  }
664  }
665  if(searchUpToQuark && !isQuarkFound) return 0;
666  if(isFromB) return 5;
667  else return 4;
668 
669 }
670 //____________________________________________________________________________
671 Double_t AliVertexingHFUtils::GetBeautyMotherPt(TClonesArray* arrayMC, AliAODMCParticle *mcPart){
673 
674  Int_t pdgGranma = 0;
675  Int_t mother = 0;
676  mother = mcPart->GetMother();
677  Int_t istep = 0;
678  Int_t abspdgGranma =0;
679  while (mother >0 ){
680  istep++;
681  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
682  if (mcGranma){
683  pdgGranma = mcGranma->GetPdgCode();
684  abspdgGranma = TMath::Abs(pdgGranma);
685  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
686  return mcGranma->Pt();
687  }
688  if(abspdgGranma==4) return -999.;
689  if(abspdgGranma==5) return -1.;
690  mother = mcGranma->GetMother();
691  }else{
692  printf("AliVertexingHFUtils::GetBeautyMotherPt: Failed casting the mother particle!");
693  break;
694  }
695  }
696  return -999.;
697 }
698 //____________________________________________________________________________
699 Int_t AliVertexingHFUtils::CheckD0Decay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
701 
702  if(label<0) return -1;
703  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
704  TParticle* part = mcEvent->Particle(label);
705  if(!part || !mcPart) return -1;
706  Int_t pdgD=part->GetPdgCode();
707  if(TMath::Abs(pdgD)!=421) return -1;
708 
709  Int_t nDau=mcPart->GetNDaughters();
710 
711  if(nDau==2){
712  Int_t daughter0 = mcPart->GetDaughterFirst();
713  Int_t daughter1 = mcPart->GetDaughterLast();
714  TParticle* mcPartDaughter0 = mcEvent->Particle(daughter0);
715  TParticle* mcPartDaughter1 = mcEvent->Particle(daughter1);
716  if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
717  arrayDauLab[0]=daughter0;
718  arrayDauLab[1]=daughter1;
719  Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
720  Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
721  if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
722  !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
723  return -1;
724  }
725  if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
726  if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
727  if((pdgCode0*pdgCode1)>0) return -1;
728  Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
729  Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
730  Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
731  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
732  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
733  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
734  return 1;
735  }
736 
737  if(nDau==3 || nDau==4){
738  Int_t nKaons=0;
739  Int_t nPions=0;
740  Double_t sumPxDau=0.;
741  Double_t sumPyDau=0.;
742  Double_t sumPzDau=0.;
743  Int_t nFoundKpi=0;
744  Int_t labelFirstDau = mcPart->GetDaughterFirst();
745  for(Int_t iDau=0; iDau<nDau; iDau++){
746  Int_t indDau = labelFirstDau+iDau;
747  if(indDau<0) return -1;
748  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
749  TParticle* dau=mcEvent->Particle(indDau);
750  if(!mcDau || !dau) return -1;
751  Int_t pdgdau=dau->GetPdgCode();
752  if(TMath::Abs(pdgdau)==321){
753  if(pdgD>0 && pdgdau>0) return -1;
754  if(pdgD<0 && pdgdau<0) return -1;
755  nKaons++;
756  sumPxDau+=dau->Px();
757  sumPyDau+=dau->Py();
758  sumPzDau+=dau->Pz();
759  arrayDauLab[nFoundKpi++]=indDau;
760  if(nFoundKpi>4) return -1;
761  }else if(TMath::Abs(pdgdau)==211){
762  nPions++;
763  sumPxDau+=dau->Px();
764  sumPyDau+=dau->Py();
765  sumPzDau+=dau->Pz();
766  arrayDauLab[nFoundKpi++]=indDau;
767  if(nFoundKpi>4) return -1;
768  }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
769  Int_t nResDau=mcDau->GetNDaughters();
770  if(nResDau!=2) return -1;
771  Int_t indFirstResDau=mcDau->GetDaughterFirst();
772  for(Int_t resDau=0; resDau<2; resDau++){
773  Int_t indResDau=indFirstResDau+resDau;
774  if(indResDau<0) return -1;
775  TParticle* resdau=mcEvent->Particle(indResDau);
776  if(!resdau) return -1;
777  Int_t pdgresdau=resdau->GetPdgCode();
778  if(TMath::Abs(pdgresdau)==321){
779  if(pdgD>0 && pdgresdau>0) return -1;
780  if(pdgD<0 && pdgresdau<0) return -1;
781  nKaons++;
782  sumPxDau+=resdau->Px();
783  sumPyDau+=resdau->Py();
784  sumPzDau+=resdau->Pz();
785  arrayDauLab[nFoundKpi++]=indResDau;
786  if(nFoundKpi>4) return -1;
787  }
788  if(TMath::Abs(pdgresdau)==211){
789  nPions++;
790  sumPxDau+=resdau->Px();
791  sumPyDau+=resdau->Py();
792  sumPzDau+=resdau->Pz();
793  arrayDauLab[nFoundKpi++]=indResDau;
794  if(nFoundKpi>4) return -1;
795  }
796  }
797  }else{
798  return -1;
799  }
800  }
801  if(nPions!=3) return -1;
802  if(nKaons!=1) return -1;
803  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -1;
804  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -1;
805  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -1;
806  return 2;
807  }
808  return -1;
809 }
810 
811 //____________________________________________________________________________
812 Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
814 
815  Int_t pdgD=mcPart->GetPdgCode();
816  if(TMath::Abs(pdgD)!=421) return -1;
817 
818  Int_t nDau=mcPart->GetNDaughters();
819 
820  if(nDau==2){
821  Int_t daughter0 = mcPart->GetDaughterLabel(0);
822  Int_t daughter1 = mcPart->GetDaughterLabel(1);
823  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
824  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
825  if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
826  arrayDauLab[0]=daughter0;
827  arrayDauLab[1]=daughter1;
828  Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
829  Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
830  if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
831  !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
832  return -1;
833  }
834  if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
835  if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
836  if((pdgCode0*pdgCode1)>0) return -1;
837  Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
838  Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
839  Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
840  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
841  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
842  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
843  return 1;
844  }
845 
846  if(nDau==3 || nDau==4){
847  Int_t nKaons=0;
848  Int_t nPions=0;
849  Double_t sumPxDau=0.;
850  Double_t sumPyDau=0.;
851  Double_t sumPzDau=0.;
852  Int_t nFoundKpi=0;
853  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
854  for(Int_t iDau=0; iDau<nDau; iDau++){
855  Int_t indDau = labelFirstDau+iDau;
856  if(indDau<0) return -1;
857  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
858  if(!dau) return -1;
859  Int_t pdgdau=dau->GetPdgCode();
860  if(TMath::Abs(pdgdau)==321){
861  if(pdgD>0 && pdgdau>0) return -1;
862  if(pdgD<0 && pdgdau<0) return -1;
863  nKaons++;
864  sumPxDau+=dau->Px();
865  sumPyDau+=dau->Py();
866  sumPzDau+=dau->Pz();
867  arrayDauLab[nFoundKpi++]=indDau;
868  if(nFoundKpi>4) return -1;
869  }else if(TMath::Abs(pdgdau)==211){
870  nPions++;
871  sumPxDau+=dau->Px();
872  sumPyDau+=dau->Py();
873  sumPzDau+=dau->Pz();
874  arrayDauLab[nFoundKpi++]=indDau;
875  if(nFoundKpi>4) return -1;
876  }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
877  Int_t nResDau=dau->GetNDaughters();
878  if(nResDau!=2) return -1;
879  Int_t indFirstResDau=dau->GetDaughterLabel(0);
880  for(Int_t resDau=0; resDau<2; resDau++){
881  Int_t indResDau=indFirstResDau+resDau;
882  if(indResDau<0) return -1;
883  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
884  if(!resdau) return -1;
885  Int_t pdgresdau=resdau->GetPdgCode();
886  if(TMath::Abs(pdgresdau)==321){
887  if(pdgD>0 && pdgresdau>0) return -1;
888  if(pdgD<0 && pdgresdau<0) return -1;
889  nKaons++;
890  sumPxDau+=resdau->Px();
891  sumPyDau+=resdau->Py();
892  sumPzDau+=resdau->Pz();
893  arrayDauLab[nFoundKpi++]=indResDau;
894  if(nFoundKpi>4) return -1;
895  }
896  if(TMath::Abs(pdgresdau)==211){
897  nPions++;
898  sumPxDau+=resdau->Px();
899  sumPyDau+=resdau->Py();
900  sumPzDau+=resdau->Pz();
901  arrayDauLab[nFoundKpi++]=indResDau;
902  if(nFoundKpi>4) return -1;
903  }
904  }
905  }else{
906  return -1;
907  }
908  }
909  if(nPions!=3) return -1;
910  if(nKaons!=1) return -1;
911  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
912  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
913  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
914  return 2;
915  }
916  return -1;
917 }
918 //____________________________________________________________________________
919 Int_t AliVertexingHFUtils::CheckDplusDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
921 
922  if(label<0) return -1;
923  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
924  TParticle* part = mcEvent->Particle(label);
925  if(!part || !mcPart) return -1;
926  Int_t pdgD=part->GetPdgCode();
927  if(TMath::Abs(pdgD)!=411) return -1;
928 
929  Int_t nDau=mcPart->GetNDaughters();
930  Int_t labelFirstDau = mcPart->GetDaughterFirst();
931  Int_t nKaons=0;
932  Int_t nPions=0;
933  Double_t sumPxDau=0.;
934  Double_t sumPyDau=0.;
935  Double_t sumPzDau=0.;
936  Int_t nFoundKpi=0;
937 
938  if(nDau==3 || nDau==2){
939  for(Int_t iDau=0; iDau<nDau; iDau++){
940  Int_t indDau = labelFirstDau+iDau;
941  if(indDau<0) return -1;
942  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
943  TParticle* dau=mcEvent->Particle(indDau);
944  if(!mcDau || !dau) return -1;
945  Int_t pdgdau=dau->GetPdgCode();
946  if(TMath::Abs(pdgdau)==321){
947  if(pdgD*pdgdau>0) return -1;
948  nKaons++;
949  sumPxDau+=dau->Px();
950  sumPyDau+=dau->Py();
951  sumPzDau+=dau->Pz();
952  arrayDauLab[nFoundKpi++]=indDau;
953  if(nFoundKpi>3) return -1;
954  }else if(TMath::Abs(pdgdau)==211){
955  if(pdgD*pdgdau<0) return -1;
956  nPions++;
957  sumPxDau+=dau->Px();
958  sumPyDau+=dau->Py();
959  sumPzDau+=dau->Pz();
960  arrayDauLab[nFoundKpi++]=indDau;
961  if(nFoundKpi>3) return -1;
962  }else if(TMath::Abs(pdgdau)==313){
963  Int_t nResDau=mcDau->GetNDaughters();
964  if(nResDau!=2) return -1;
965  Int_t indFirstResDau=mcDau->GetDaughterFirst();
966  for(Int_t resDau=0; resDau<2; resDau++){
967  Int_t indResDau=indFirstResDau+resDau;
968  if(indResDau<0) return -1;
969  TParticle* resdau=mcEvent->Particle(indResDau);
970  if(!resdau) return -1;
971  Int_t pdgresdau=resdau->GetPdgCode();
972  if(TMath::Abs(pdgresdau)==321){
973  if(pdgD*pdgresdau>0) return -1;
974  sumPxDau+=resdau->Px();
975  sumPyDau+=resdau->Py();
976  sumPzDau+=resdau->Pz();
977  nKaons++;
978  arrayDauLab[nFoundKpi++]=indResDau;
979  if(nFoundKpi>3) return -1;
980  }
981  if(TMath::Abs(pdgresdau)==211){
982  if(pdgD*pdgresdau<0) return -1;
983  sumPxDau+=resdau->Px();
984  sumPyDau+=resdau->Py();
985  sumPzDau+=resdau->Pz();
986  nPions++;
987  arrayDauLab[nFoundKpi++]=indResDau;
988  if(nFoundKpi>3) return -1;
989  }
990  }
991  }else{
992  return -1;
993  }
994  }
995  if(nPions!=2) return -1;
996  if(nKaons!=1) return -1;
997  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
998  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
999  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1000  if(nDau==3) return 1;
1001  else if(nDau==2) return 2;
1002  }
1003 
1004  return -1;
1005 
1006 }
1007 
1008 //____________________________________________________________________________
1009 Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1011 
1012  Int_t pdgD=mcPart->GetPdgCode();
1013  if(TMath::Abs(pdgD)!=411) return -1;
1014 
1015  Int_t nDau=mcPart->GetNDaughters();
1016  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1017  Int_t nKaons=0;
1018  Int_t nPions=0;
1019  Double_t sumPxDau=0.;
1020  Double_t sumPyDau=0.;
1021  Double_t sumPzDau=0.;
1022  Int_t nFoundKpi=0;
1023 
1024  if(nDau==3 || nDau==2){
1025  for(Int_t iDau=0; iDau<nDau; iDau++){
1026  Int_t indDau = labelFirstDau+iDau;
1027  if(indDau<0) return -1;
1028  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1029  if(!dau) return -1;
1030  Int_t pdgdau=dau->GetPdgCode();
1031  if(TMath::Abs(pdgdau)==321){
1032  if(pdgD*pdgdau>0) return -1;
1033  nKaons++;
1034  sumPxDau+=dau->Px();
1035  sumPyDau+=dau->Py();
1036  sumPzDau+=dau->Pz();
1037  arrayDauLab[nFoundKpi++]=indDau;
1038  if(nFoundKpi>3) return -1;
1039  }else if(TMath::Abs(pdgdau)==211){
1040  if(pdgD*pdgdau<0) return -1;
1041  nPions++;
1042  sumPxDau+=dau->Px();
1043  sumPyDau+=dau->Py();
1044  sumPzDau+=dau->Pz();
1045  arrayDauLab[nFoundKpi++]=indDau;
1046  if(nFoundKpi>3) return -1;
1047  }else if(TMath::Abs(pdgdau)==313){
1048  Int_t nResDau=dau->GetNDaughters();
1049  if(nResDau!=2) return -1;
1050  Int_t indFirstResDau=dau->GetDaughterLabel(0);
1051  for(Int_t resDau=0; resDau<2; resDau++){
1052  Int_t indResDau=indFirstResDau+resDau;
1053  if(indResDau<0) return -1;
1054  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1055  if(!resdau) return -1;
1056  Int_t pdgresdau=resdau->GetPdgCode();
1057  if(TMath::Abs(pdgresdau)==321){
1058  if(pdgD*pdgresdau>0) return -1;
1059  sumPxDau+=resdau->Px();
1060  sumPyDau+=resdau->Py();
1061  sumPzDau+=resdau->Pz();
1062  nKaons++;
1063  arrayDauLab[nFoundKpi++]=indResDau;
1064  if(nFoundKpi>3) return -1;
1065  }
1066  if(TMath::Abs(pdgresdau)==211){
1067  if(pdgD*pdgresdau<0) return -1;
1068  sumPxDau+=resdau->Px();
1069  sumPyDau+=resdau->Py();
1070  sumPzDau+=resdau->Pz();
1071  nPions++;
1072  arrayDauLab[nFoundKpi++]=indResDau;
1073  if(nFoundKpi>3) return -1;
1074  }
1075  }
1076  }else{
1077  return -1;
1078  }
1079  }
1080  if(nPions!=2) return -1;
1081  if(nKaons!=1) return -1;
1082  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1083  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1084  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1085  if(nDau==3) return 1;
1086  else if(nDau==2) return 2;
1087  }
1088 
1089  return -1;
1090 
1091 }
1092 //____________________________________________________________________________
1093 Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1095 
1096  if(label<0) return -1;
1097  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1098  TParticle* part = mcEvent->Particle(label);
1099  if(!part || !mcPart) return -1;
1100  Int_t pdgD=part->GetPdgCode();
1101  if(TMath::Abs(pdgD)!=411) return -1;
1102 
1103  Int_t nDau=mcPart->GetNDaughters();
1104  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1105  Int_t nKaons=0;
1106  Int_t nPions=0;
1107  Double_t sumPxDau=0.;
1108  Double_t sumPyDau=0.;
1109  Double_t sumPzDau=0.;
1110  Int_t nFoundKpi=0;
1111  Bool_t isPhi=kFALSE;
1112  Bool_t isk0st=kFALSE;
1113 
1114  if(nDau==3 || nDau==2){
1115  for(Int_t iDau=0; iDau<nDau; iDau++){
1116  Int_t indDau = labelFirstDau+iDau;
1117  if(indDau<0) return -1;
1118  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1119  TParticle* dau=mcEvent->Particle(indDau);
1120  if(!mcDau || !dau) return -1;
1121  Int_t pdgdau=dau->GetPdgCode();
1122  if(TMath::Abs(pdgdau)==321){
1123  nKaons++;
1124  sumPxDau+=dau->Px();
1125  sumPyDau+=dau->Py();
1126  sumPzDau+=dau->Pz();
1127  arrayDauLab[nFoundKpi++]=indDau;
1128  if(nFoundKpi>3) return -1;
1129  }else if(TMath::Abs(pdgdau)==211){
1130  nPions++;
1131  sumPxDau+=dau->Px();
1132  sumPyDau+=dau->Py();
1133  sumPzDau+=dau->Pz();
1134  arrayDauLab[nFoundKpi++]=indDau;
1135  if(nFoundKpi>3) return -1;
1136  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1137  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1138  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1139  Int_t nResDau=mcDau->GetNDaughters();
1140  if(nResDau!=2) return -1;
1141  Int_t indFirstResDau=mcDau->GetDaughterFirst();
1142  for(Int_t resDau=0; resDau<2; resDau++){
1143  Int_t indResDau=indFirstResDau+resDau;
1144  if(indResDau<0) return -1;
1145  TParticle* resdau=mcEvent->Particle(indResDau);
1146  if(!resdau) return -1;
1147  Int_t pdgresdau=resdau->GetPdgCode();
1148  if(TMath::Abs(pdgresdau)==321){
1149  sumPxDau+=resdau->Px();
1150  sumPyDau+=resdau->Py();
1151  sumPzDau+=resdau->Pz();
1152  nKaons++;
1153  arrayDauLab[nFoundKpi++]=indResDau;
1154  if(nFoundKpi>3) return -1;
1155  }
1156  if(TMath::Abs(pdgresdau)==211){
1157  sumPxDau+=resdau->Px();
1158  sumPyDau+=resdau->Py();
1159  sumPzDau+=resdau->Pz();
1160  nPions++;
1161  arrayDauLab[nFoundKpi++]=indResDau;
1162  if(nFoundKpi>3) return -1;
1163  }
1164  }
1165  }else{
1166  return -1;
1167  }
1168  }
1169  if(nPions!=1) return -1;
1170  if(nKaons!=2) return -1;
1171  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1172  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1173  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1174  if(nDau==3) return 3;
1175  else if(nDau==2){
1176  if(isk0st) return 2;
1177  if(isPhi) return 1;
1178  }
1179  }
1180 
1181  return -1;
1182 
1183 }
1184 //____________________________________________________________________________
1185 Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1187 
1188  Int_t pdgD=mcPart->GetPdgCode();
1189  if(TMath::Abs(pdgD)!=411) return -1;
1190 
1191  Int_t nDau=mcPart->GetNDaughters();
1192  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1193  Int_t nKaons=0;
1194  Int_t nPions=0;
1195  Double_t sumPxDau=0.;
1196  Double_t sumPyDau=0.;
1197  Double_t sumPzDau=0.;
1198  Int_t nFoundKpi=0;
1199  Bool_t isPhi=kFALSE;
1200  Bool_t isk0st=kFALSE;
1201 
1202  if(nDau==3 || nDau==2){
1203  for(Int_t iDau=0; iDau<nDau; iDau++){
1204  Int_t indDau = labelFirstDau+iDau;
1205  if(indDau<0) return -1;
1206  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1207  if(!dau) return -1;
1208  Int_t pdgdau=dau->GetPdgCode();
1209  if(TMath::Abs(pdgdau)==321){
1210  nKaons++;
1211  sumPxDau+=dau->Px();
1212  sumPyDau+=dau->Py();
1213  sumPzDau+=dau->Pz();
1214  arrayDauLab[nFoundKpi++]=indDau;
1215  if(nFoundKpi>3) return -1;
1216  }else if(TMath::Abs(pdgdau)==211){
1217  nPions++;
1218  sumPxDau+=dau->Px();
1219  sumPyDau+=dau->Py();
1220  sumPzDau+=dau->Pz();
1221  arrayDauLab[nFoundKpi++]=indDau;
1222  if(nFoundKpi>3) return -1;
1223  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1224  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1225  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1226  Int_t nResDau=dau->GetNDaughters();
1227  if(nResDau!=2) return -1;
1228  Int_t indFirstResDau=dau->GetDaughterLabel(0);
1229  for(Int_t resDau=0; resDau<2; resDau++){
1230  Int_t indResDau=indFirstResDau+resDau;
1231  if(indResDau<0) return -1;
1232  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1233  if(!resdau) return -1;
1234  Int_t pdgresdau=resdau->GetPdgCode();
1235  if(TMath::Abs(pdgresdau)==321){
1236  sumPxDau+=resdau->Px();
1237  sumPyDau+=resdau->Py();
1238  sumPzDau+=resdau->Pz();
1239  nKaons++;
1240  arrayDauLab[nFoundKpi++]=indResDau;
1241  if(nFoundKpi>3) return -1;
1242  }
1243  if(TMath::Abs(pdgresdau)==211){
1244  sumPxDau+=resdau->Px();
1245  sumPyDau+=resdau->Py();
1246  sumPzDau+=resdau->Pz();
1247  nPions++;
1248  arrayDauLab[nFoundKpi++]=indResDau;
1249  if(nFoundKpi>3) return -1;
1250  }
1251  }
1252  }else{
1253  return -1;
1254  }
1255  }
1256  if(nPions!=1) return -1;
1257  if(nKaons!=2) return -1;
1258  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1259  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1260  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1261  if(nDau==3) return 3;
1262  else if(nDau==2){
1263  if(isk0st) return 2;
1264  if(isPhi) return 1;
1265  }
1266  }
1267 
1268  return -1;
1269 
1270 }
1271 //____________________________________________________________________________
1272 Int_t AliVertexingHFUtils::CheckDplusK0spiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1274 
1275  if(label<0) return -1;
1276  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1277  TParticle* part = mcEvent->Particle(label);
1278  if(!part || !mcPart) return -1;
1279  Int_t pdgD=part->GetPdgCode();
1280  if(TMath::Abs(pdgD)!=411) return -1;
1281 
1282  Int_t nDau=mcPart->GetNDaughters();
1283  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1284  Int_t nPions=0;
1285  Double_t sumPxDau=0.;
1286  Double_t sumPyDau=0.;
1287  Double_t sumPzDau=0.;
1288  Int_t nFoundpi=0;
1289 
1290  Int_t codeV0=-1;
1291  if(nDau==2){
1292  for(Int_t iDau=0; iDau<nDau; iDau++){
1293  Int_t indDau = labelFirstDau+iDau;
1294  if(indDau<0) return -1;
1295  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1296  TParticle* dau=mcEvent->Particle(indDau);
1297  if(!mcDau || !dau) return -1;
1298  Int_t pdgdau=dau->GetPdgCode();
1299  if(TMath::Abs(pdgdau)==211){
1300  nPions++;
1301  sumPxDau+=dau->Px();
1302  sumPyDau+=dau->Py();
1303  sumPzDau+=dau->Pz();
1304  arrayDauLab[nFoundpi++]=indDau;
1305  if(nFoundpi>3) return -1;
1306  }else if(TMath::Abs(pdgdau)==311){
1307  codeV0=TMath::Abs(pdgdau);
1308  TParticle* v0=dau;
1309  AliMCParticle* mcV0=mcDau;
1310  if(codeV0==311){
1311  Int_t nK0Dau=mcDau->GetNDaughters();
1312  if(nK0Dau!=1) return -1;
1313  Int_t indK0s=mcDau->GetDaughterFirst();
1314  if(indK0s<0) return -1;
1315  v0=mcEvent->Particle(indK0s);
1316  mcV0=(AliMCParticle*)mcEvent->GetTrack(indK0s);
1317  if(!v0 || !mcV0) return -1;
1318  Int_t pdgK0sl=v0->GetPdgCode();
1319  codeV0=TMath::Abs(pdgK0sl);
1320  }
1321  Int_t nV0Dau=mcV0->GetNDaughters();
1322  if(nV0Dau!=2) return -1;
1323  Int_t indFirstV0Dau=mcV0->GetDaughterFirst();
1324  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1325  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1326  if(indV0Dau<0) return -1;
1327  TParticle* v0dau=mcEvent->Particle(indV0Dau);
1328  if(!v0dau) return -1;
1329  Int_t pdgv0dau=v0dau->GetPdgCode();
1330  if(TMath::Abs(pdgv0dau)==211){
1331  sumPxDau+=v0dau->Px();
1332  sumPyDau+=v0dau->Py();
1333  sumPzDau+=v0dau->Pz();
1334  nPions++;
1335  arrayDauLab[nFoundpi++]=indV0Dau;
1336  if(nFoundpi>3) return -1;
1337  }
1338  }
1339  }else{
1340  return -1;
1341  }
1342  }
1343  if(nPions!=3) return -1;
1344  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1345  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1346  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1347  if(codeV0==310) return 1;
1348  }
1349  return -1;
1350 
1351 }
1352 
1353 
1354 //____________________________________________________________________________
1355 Int_t AliVertexingHFUtils::CheckDsDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1357 
1358  if(label<0) return -1;
1359  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1360  TParticle* part = mcEvent->Particle(label);
1361  if(!part || !mcPart) return -1;
1362  Int_t pdgD=part->GetPdgCode();
1363  if(TMath::Abs(pdgD)!=431) return -1;
1364 
1365  Int_t nDau=mcPart->GetNDaughters();
1366  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1367  Int_t nKaons=0;
1368  Int_t nPions=0;
1369  Double_t sumPxDau=0.;
1370  Double_t sumPyDau=0.;
1371  Double_t sumPzDau=0.;
1372  Int_t nFoundKpi=0;
1373  Bool_t isPhi=kFALSE;
1374  Bool_t isk0st=kFALSE;
1375  Bool_t isf0=kFALSE;
1376 
1377  if(nDau==3 || nDau==2){
1378  for(Int_t iDau=0; iDau<nDau; iDau++){
1379  Int_t indDau = labelFirstDau+iDau;
1380  if(indDau<0) return -1;
1381  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1382  TParticle* dau=mcEvent->Particle(indDau);
1383  if(!mcDau || !dau) return -1;
1384  Int_t pdgdau=dau->GetPdgCode();
1385  if(TMath::Abs(pdgdau)==321){
1386  nKaons++;
1387  sumPxDau+=dau->Px();
1388  sumPyDau+=dau->Py();
1389  sumPzDau+=dau->Pz();
1390  arrayDauLab[nFoundKpi++]=indDau;
1391  if(nFoundKpi>3) return -1;
1392  }else if(TMath::Abs(pdgdau)==211){
1393  nPions++;
1394  sumPxDau+=dau->Px();
1395  sumPyDau+=dau->Py();
1396  sumPzDau+=dau->Pz();
1397  arrayDauLab[nFoundKpi++]=indDau;
1398  if(nFoundKpi>3) return -1;
1399  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333 || TMath::Abs(pdgdau)==9010221){
1400  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1401  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1402  if(TMath::Abs(pdgdau)==9010221) isf0=kTRUE;
1403  Int_t nResDau=mcDau->GetNDaughters();
1404  if(nResDau!=2) return -1;
1405  Int_t indFirstResDau=mcDau->GetDaughterFirst();
1406  for(Int_t resDau=0; resDau<2; resDau++){
1407  Int_t indResDau=indFirstResDau+resDau;
1408  if(indResDau<0) return -1;
1409  TParticle* resdau=mcEvent->Particle(indResDau);
1410  if(!resdau) return -1;
1411  Int_t pdgresdau=resdau->GetPdgCode();
1412  if(TMath::Abs(pdgresdau)==321){
1413  sumPxDau+=resdau->Px();
1414  sumPyDau+=resdau->Py();
1415  sumPzDau+=resdau->Pz();
1416  nKaons++;
1417  arrayDauLab[nFoundKpi++]=indResDau;
1418  if(nFoundKpi>3) return -1;
1419  }
1420  if(TMath::Abs(pdgresdau)==211){
1421  sumPxDau+=resdau->Px();
1422  sumPyDau+=resdau->Py();
1423  sumPzDau+=resdau->Pz();
1424  nPions++;
1425  arrayDauLab[nFoundKpi++]=indResDau;
1426  if(nFoundKpi>3) return -1;
1427  }
1428  }
1429  }else{
1430  return -1;
1431  }
1432  }
1433  if(nPions!=1) return -1;
1434  if(nKaons!=2) return -1;
1435  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1436  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1437  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1438  if(nDau==3) return 3;
1439  else if(nDau==2){
1440  if(isk0st) return 2;
1441  if(isPhi) return 1;
1442  if(isf0) return 4;
1443  }
1444  }
1445 
1446  return -1;
1447 }
1448 //____________________________________________________________________________
1449 Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1451 
1452  Int_t pdgD=mcPart->GetPdgCode();
1453  if(TMath::Abs(pdgD)!=431) return -1;
1454 
1455  Int_t nDau=mcPart->GetNDaughters();
1456  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1457  Int_t nKaons=0;
1458  Int_t nPions=0;
1459  Double_t sumPxDau=0.;
1460  Double_t sumPyDau=0.;
1461  Double_t sumPzDau=0.;
1462  Int_t nFoundKpi=0;
1463  Bool_t isPhi=kFALSE;
1464  Bool_t isk0st=kFALSE;
1465  Bool_t isf0=kFALSE;
1466 
1467  if(nDau==3 || nDau==2){
1468  for(Int_t iDau=0; iDau<nDau; iDau++){
1469  Int_t indDau = labelFirstDau+iDau;
1470  if(indDau<0) return -1;
1471  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1472  if(!dau) return -1;
1473  Int_t pdgdau=dau->GetPdgCode();
1474  if(TMath::Abs(pdgdau)==321){
1475  nKaons++;
1476  sumPxDau+=dau->Px();
1477  sumPyDau+=dau->Py();
1478  sumPzDau+=dau->Pz();
1479  arrayDauLab[nFoundKpi++]=indDau;
1480  if(nFoundKpi>3) return -1;
1481  }else if(TMath::Abs(pdgdau)==211){
1482  nPions++;
1483  sumPxDau+=dau->Px();
1484  sumPyDau+=dau->Py();
1485  sumPzDau+=dau->Pz();
1486  arrayDauLab[nFoundKpi++]=indDau;
1487  if(nFoundKpi>3) return -1;
1488  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333 || TMath::Abs(pdgdau)==9010221){
1489  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1490  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1491  if(TMath::Abs(pdgdau)==9010221) isf0=kTRUE;
1492  Int_t nResDau=dau->GetNDaughters();
1493  if(nResDau!=2) return -1;
1494  Int_t indFirstResDau=dau->GetDaughterLabel(0);
1495  for(Int_t resDau=0; resDau<2; resDau++){
1496  Int_t indResDau=indFirstResDau+resDau;
1497  if(indResDau<0) return -1;
1498  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1499  if(!resdau) return -1;
1500  Int_t pdgresdau=resdau->GetPdgCode();
1501  if(TMath::Abs(pdgresdau)==321){
1502  sumPxDau+=resdau->Px();
1503  sumPyDau+=resdau->Py();
1504  sumPzDau+=resdau->Pz();
1505  nKaons++;
1506  arrayDauLab[nFoundKpi++]=indResDau;
1507  if(nFoundKpi>3) return -1;
1508  }
1509  if(TMath::Abs(pdgresdau)==211){
1510  sumPxDau+=resdau->Px();
1511  sumPyDau+=resdau->Py();
1512  sumPzDau+=resdau->Pz();
1513  nPions++;
1514  arrayDauLab[nFoundKpi++]=indResDau;
1515  if(nFoundKpi>3) return -1;
1516  }
1517  }
1518  }else{
1519  return -1;
1520  }
1521  }
1522  if(nPions!=1) return -1;
1523  if(nKaons!=2) return -1;
1524  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1525  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1526  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1527  if(nDau==3) return 3;
1528  else if(nDau==2){
1529  if(isk0st) return 2;
1530  if(isPhi) return 1;
1531  if(isf0) return 4;
1532  }
1533  }
1534 
1535  return -1;
1536 }
1537 //____________________________________________________________________________
1538 Int_t AliVertexingHFUtils::CheckDsK0sKDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1540 
1541  if(label<0) return -1;
1542  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1543  TParticle* part = mcEvent->Particle(label);
1544  if(!part || !mcPart) return -1;
1545  Int_t pdgD=part->GetPdgCode();
1546  if(TMath::Abs(pdgD)!=431) return -1;
1547 
1548  Int_t nDau=mcPart->GetNDaughters();
1549  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1550  Int_t nPions=0;
1551  Int_t nKaons=0;
1552  Double_t sumPxDau=0.;
1553  Double_t sumPyDau=0.;
1554  Double_t sumPzDau=0.;
1555  Int_t nFoundKpi=0;
1556 
1557  Int_t codeV0=-1;
1558  if(nDau==2){
1559  for(Int_t iDau=0; iDau<nDau; iDau++){
1560  Int_t indDau = labelFirstDau+iDau;
1561  if(indDau<0) return -1;
1562  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1563  TParticle* dau=mcEvent->Particle(indDau);
1564  if(!mcDau || !dau) return -1;
1565  Int_t pdgdau=dau->GetPdgCode();
1566  if(TMath::Abs(pdgdau)==211){
1567  nPions++;
1568  sumPxDau+=dau->Px();
1569  sumPyDau+=dau->Py();
1570  sumPzDau+=dau->Pz();
1571  arrayDauLab[nFoundKpi++]=indDau;
1572  if(nFoundKpi>3) return -1;
1573  }else if(TMath::Abs(pdgdau)==321){
1574  nKaons++;
1575  sumPxDau+=dau->Px();
1576  sumPyDau+=dau->Py();
1577  sumPzDau+=dau->Pz();
1578  arrayDauLab[nFoundKpi++]=indDau;
1579  if(nFoundKpi>3) return -1;
1580  }else if(TMath::Abs(pdgdau)==311){
1581  codeV0=TMath::Abs(pdgdau);
1582  TParticle* v0=dau;
1583  AliMCParticle* mcV0=mcDau;
1584  if(codeV0==311){
1585  Int_t nK0Dau=mcDau->GetNDaughters();
1586  if(nK0Dau!=1) return -1;
1587  Int_t indK0s=mcDau->GetDaughterFirst();
1588  if(indK0s<0) return -1;
1589  v0=mcEvent->Particle(indK0s);
1590  mcV0=(AliMCParticle*)mcEvent->GetTrack(indK0s);
1591  if(!v0 || !mcV0) return -1;
1592  Int_t pdgK0sl=v0->GetPdgCode();
1593  codeV0=TMath::Abs(pdgK0sl);
1594  }
1595  Int_t nV0Dau=mcV0->GetNDaughters();
1596  if(nV0Dau!=2) return -1;
1597  Int_t indFirstV0Dau=mcV0->GetDaughterFirst();
1598  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1599  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1600  if(indV0Dau<0) return -1;
1601  TParticle* v0dau=mcEvent->Particle(indV0Dau);
1602  if(!v0dau) return -1;
1603  Int_t pdgv0dau=v0dau->GetPdgCode();
1604  if(TMath::Abs(pdgv0dau)==211){
1605  sumPxDau+=v0dau->Px();
1606  sumPyDau+=v0dau->Py();
1607  sumPzDau+=v0dau->Pz();
1608  nPions++;
1609  arrayDauLab[nFoundKpi++]=indV0Dau;
1610  if(nFoundKpi>3) return -1;
1611  }
1612  }
1613  }else{
1614  return -1;
1615  }
1616  }
1617  if(nPions!=2) return -1;
1618  if(nKaons!=1) return -1;
1619  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1620  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1621  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1622  if(codeV0==310) return 1;
1623  }
1624  return -1;
1625 
1626 }
1627 //____________________________________________________________________________
1628 Int_t AliVertexingHFUtils::CheckDstarDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1630 
1631  if(label<0) return -1;
1632  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1633  TParticle* part = mcEvent->Particle(label);
1634  if(!part || !mcPart) return -1;
1635  Int_t pdgD=part->GetPdgCode();
1636  if(TMath::Abs(pdgD)!=413) return -1;
1637 
1638  Int_t nDau=mcPart->GetNDaughters();
1639  if(nDau!=2) return -1;
1640 
1641  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1642  Int_t nKaons=0;
1643  Int_t nPions=0;
1644  Double_t sumPxDau=0.;
1645  Double_t sumPyDau=0.;
1646  Double_t sumPzDau=0.;
1647  Int_t nFoundKpi=0;
1648 
1649  for(Int_t iDau=0; iDau<nDau; iDau++){
1650  Int_t indDau = labelFirstDau+iDau;
1651  if(indDau<0) return -1;
1652  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1653  TParticle* dau=mcEvent->Particle(indDau);
1654  if(!mcDau || !dau) return -1;
1655  Int_t pdgdau=dau->GetPdgCode();
1656  if(TMath::Abs(pdgdau)==421){
1657  Int_t nResDau=mcDau->GetNDaughters();
1658  if(nResDau!=2) return -1;
1659  Int_t indFirstResDau=mcDau->GetDaughterFirst();
1660  for(Int_t resDau=0; resDau<2; resDau++){
1661  Int_t indResDau=indFirstResDau+resDau;
1662  if(indResDau<0) return -1;
1663  TParticle* resdau=mcEvent->Particle(indResDau);
1664  if(!resdau) return -1;
1665  Int_t pdgresdau=resdau->GetPdgCode();
1666  if(TMath::Abs(pdgresdau)==321){
1667  if(pdgD*pdgresdau>0) return -1;
1668  sumPxDau+=resdau->Px();
1669  sumPyDau+=resdau->Py();
1670  sumPzDau+=resdau->Pz();
1671  nKaons++;
1672  arrayDauLab[nFoundKpi++]=indResDau;
1673  if(nFoundKpi>3) return -1;
1674  }
1675  if(TMath::Abs(pdgresdau)==211){
1676  if(pdgD*pdgresdau<0) return -1;
1677  sumPxDau+=resdau->Px();
1678  sumPyDau+=resdau->Py();
1679  sumPzDau+=resdau->Pz();
1680  nPions++;
1681  arrayDauLab[nFoundKpi++]=indResDau;
1682  if(nFoundKpi>3) return -1;
1683  }
1684  }
1685  }else if(TMath::Abs(pdgdau)==211){
1686  if(pdgD*pdgdau<0) return -1;
1687  nPions++;
1688  sumPxDau+=dau->Px();
1689  sumPyDau+=dau->Py();
1690  sumPzDau+=dau->Pz();
1691  arrayDauLab[nFoundKpi++]=indDau;
1692  if(nFoundKpi>3) return -1;
1693  }
1694  }
1695 
1696  if(nPions!=2) return -1;
1697  if(nKaons!=1) return -1;
1698  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1699  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1700  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1701  return 1;
1702 
1703 }
1704 
1705 //____________________________________________________________________________
1706 Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1708 
1709  Int_t pdgD=mcPart->GetPdgCode();
1710  if(TMath::Abs(pdgD)!=413) return -1;
1711 
1712  Int_t nDau=mcPart->GetNDaughters();
1713  if(nDau!=2) return -1;
1714 
1715  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1716  Int_t nKaons=0;
1717  Int_t nPions=0;
1718  Double_t sumPxDau=0.;
1719  Double_t sumPyDau=0.;
1720  Double_t sumPzDau=0.;
1721  Int_t nFoundKpi=0;
1722 
1723  for(Int_t iDau=0; iDau<nDau; iDau++){
1724  Int_t indDau = labelFirstDau+iDau;
1725  if(indDau<0) return -1;
1726  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1727  if(!dau) return -1;
1728  Int_t pdgdau=dau->GetPdgCode();
1729  if(TMath::Abs(pdgdau)==421){
1730  Int_t nResDau=dau->GetNDaughters();
1731  if(nResDau!=2) return -1;
1732  Int_t indFirstResDau=dau->GetDaughterLabel(0);
1733  for(Int_t resDau=0; resDau<2; resDau++){
1734  Int_t indResDau=indFirstResDau+resDau;
1735  if(indResDau<0) return -1;
1736  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1737  if(!resdau) return -1;
1738  Int_t pdgresdau=resdau->GetPdgCode();
1739  if(TMath::Abs(pdgresdau)==321){
1740  if(pdgD*pdgresdau>0) return -1;
1741  sumPxDau+=resdau->Px();
1742  sumPyDau+=resdau->Py();
1743  sumPzDau+=resdau->Pz();
1744  nKaons++;
1745  arrayDauLab[nFoundKpi++]=indResDau;
1746  if(nFoundKpi>3) return -1;
1747  }
1748  if(TMath::Abs(pdgresdau)==211){
1749  if(pdgD*pdgresdau<0) return -1;
1750  sumPxDau+=resdau->Px();
1751  sumPyDau+=resdau->Py();
1752  sumPzDau+=resdau->Pz();
1753  nPions++;
1754  arrayDauLab[nFoundKpi++]=indResDau;
1755  if(nFoundKpi>3) return -1;
1756  }
1757  }
1758  }else if(TMath::Abs(pdgdau)==211){
1759  if(pdgD*pdgdau<0) return -1;
1760  nPions++;
1761  sumPxDau+=dau->Px();
1762  sumPyDau+=dau->Py();
1763  sumPzDau+=dau->Pz();
1764  arrayDauLab[nFoundKpi++]=indDau;
1765  if(nFoundKpi>3) return -1;
1766  }
1767  }
1768 
1769  if(nPions!=2) return -1;
1770  if(nKaons!=1) return -1;
1771  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1772  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1773  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1774  return 1;
1775 
1776 }
1777 
1778 //____________________________________________________________________________
1779 Int_t AliVertexingHFUtils::CheckLcpKpiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1781 
1782  if(label<0) return -1;
1783  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1784  TParticle* part = mcEvent->Particle(label);
1785  if(!part || !mcPart) return -1;
1786  Int_t pdgD=part->GetPdgCode();
1787  if(TMath::Abs(pdgD)!=4122) return -1;
1788 
1789  Int_t nDau=mcPart->GetNDaughters();
1790  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1791  Int_t nKaons=0;
1792  Int_t nPions=0;
1793  Int_t nProtons=0;
1794  Double_t sumPxDau=0.;
1795  Double_t sumPyDau=0.;
1796  Double_t sumPzDau=0.;
1797  Int_t nFoundpKpi=0;
1798 
1799  Int_t codeRes=-1;
1800  if(nDau==3 || nDau==2){
1801  for(Int_t iDau=0; iDau<nDau; iDau++){
1802  Int_t indDau = labelFirstDau+iDau;
1803  if(indDau<0) return -1;
1804  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1805  TParticle* dau=mcEvent->Particle(indDau);
1806  if(!mcDau || !dau) return -1;
1807  Int_t pdgdau=dau->GetPdgCode();
1808  if(TMath::Abs(pdgdau)==321){
1809  nKaons++;
1810  sumPxDau+=dau->Px();
1811  sumPyDau+=dau->Py();
1812  sumPzDau+=dau->Pz();
1813  arrayDauLab[nFoundpKpi++]=indDau;
1814  if(nFoundpKpi>3) return -1;
1815  }else if(TMath::Abs(pdgdau)==211){
1816  nPions++;
1817  sumPxDau+=dau->Px();
1818  sumPyDau+=dau->Py();
1819  sumPzDau+=dau->Pz();
1820  arrayDauLab[nFoundpKpi++]=indDau;
1821  if(nFoundpKpi>3) return -1;
1822  }else if(TMath::Abs(pdgdau)==2212){
1823  nProtons++;
1824  sumPxDau+=dau->Px();
1825  sumPyDau+=dau->Py();
1826  sumPzDau+=dau->Pz();
1827  arrayDauLab[nFoundpKpi++]=indDau;
1828  if(nFoundpKpi>3) return -1;
1829  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 ||
1830  TMath::Abs(pdgdau)==2224){
1831  codeRes=TMath::Abs(pdgdau);
1832  Int_t nResDau=mcDau->GetNDaughters();
1833  if(nResDau!=2) return -1;
1834  Int_t indFirstResDau=mcDau->GetDaughterFirst();
1835  for(Int_t resDau=0; resDau<2; resDau++){
1836  Int_t indResDau=indFirstResDau+resDau;
1837  if(indResDau<0) return -1;
1838  TParticle* resdau=mcEvent->Particle(indResDau);
1839  if(!resdau) return -1;
1840  Int_t pdgresdau=resdau->GetPdgCode();
1841  if(TMath::Abs(pdgresdau)==321){
1842  sumPxDau+=resdau->Px();
1843  sumPyDau+=resdau->Py();
1844  sumPzDau+=resdau->Pz();
1845  nKaons++;
1846  arrayDauLab[nFoundpKpi++]=indResDau;
1847  if(nFoundpKpi>3) return -1;
1848  }else if(TMath::Abs(pdgresdau)==211){
1849  sumPxDau+=resdau->Px();
1850  sumPyDau+=resdau->Py();
1851  sumPzDau+=resdau->Pz();
1852  nPions++;
1853  arrayDauLab[nFoundpKpi++]=indResDau;
1854  if(nFoundpKpi>3) return -1;
1855  }else if(TMath::Abs(pdgresdau)==2212){
1856  sumPxDau+=resdau->Px();
1857  sumPyDau+=resdau->Py();
1858  sumPzDau+=resdau->Pz();
1859  nProtons++;
1860  arrayDauLab[nFoundpKpi++]=indResDau;
1861  if(nFoundpKpi>3) return -1;
1862  }
1863  }
1864  }else{
1865  return -1;
1866  }
1867  }
1868  if(nPions!=1) return -1;
1869  if(nKaons!=1) return -1;
1870  if(nProtons!=1) return -1;
1871  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1872  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1873  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1874  if(nDau==3) return 1;
1875  else if(nDau==2){
1876  if(codeRes==313) return 2;
1877  else if(codeRes==2224) return 3;
1878  else if(codeRes==3124) return 4;
1879  }
1880  }
1881  return -1;
1882 
1883 }
1884 
1885 //____________________________________________________________________________
1886 Int_t AliVertexingHFUtils::CheckLcpKpiDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1888 
1889  Int_t pdgD=mcPart->GetPdgCode();
1890  if(TMath::Abs(pdgD)!=4122) return -1;
1891 
1892  Int_t nDau=mcPart->GetNDaughters();
1893  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
1894  Int_t nKaons=0;
1895  Int_t nPions=0;
1896  Int_t nProtons=0;
1897  Double_t sumPxDau=0.;
1898  Double_t sumPyDau=0.;
1899  Double_t sumPzDau=0.;
1900  Int_t nFoundpKpi=0;
1901 
1902  Int_t codeRes=-1;
1903  if(nDau==3 || nDau==2){
1904  for(Int_t iDau=0; iDau<nDau; iDau++){
1905  Int_t indDau = labelFirstDau+iDau;
1906  if(indDau<0) return -1;
1907  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1908  if(!dau) return -1;
1909  Int_t pdgdau=dau->GetPdgCode();
1910  if(TMath::Abs(pdgdau)==321){
1911  nKaons++;
1912  sumPxDau+=dau->Px();
1913  sumPyDau+=dau->Py();
1914  sumPzDau+=dau->Pz();
1915  arrayDauLab[nFoundpKpi++]=indDau;
1916  if(nFoundpKpi>3) return -1;
1917  }else if(TMath::Abs(pdgdau)==211){
1918  nPions++;
1919  sumPxDau+=dau->Px();
1920  sumPyDau+=dau->Py();
1921  sumPzDau+=dau->Pz();
1922  arrayDauLab[nFoundpKpi++]=indDau;
1923  if(nFoundpKpi>3) return -1;
1924  }else if(TMath::Abs(pdgdau)==2212){
1925  nProtons++;
1926  sumPxDau+=dau->Px();
1927  sumPyDau+=dau->Py();
1928  sumPzDau+=dau->Pz();
1929  arrayDauLab[nFoundpKpi++]=indDau;
1930  if(nFoundpKpi>3) return -1;
1931  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 ||
1932  TMath::Abs(pdgdau)==2224){
1933  codeRes=TMath::Abs(pdgdau);
1934  Int_t nResDau=dau->GetNDaughters();
1935  if(nResDau!=2) return -1;
1936  Int_t indFirstResDau=dau->GetDaughterLabel(0);
1937  for(Int_t resDau=0; resDau<2; resDau++){
1938  Int_t indResDau=indFirstResDau+resDau;
1939  if(indResDau<0) return -1;
1940  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1941  if(!resdau) return -1;
1942  Int_t pdgresdau=resdau->GetPdgCode();
1943  if(TMath::Abs(pdgresdau)==321){
1944  sumPxDau+=resdau->Px();
1945  sumPyDau+=resdau->Py();
1946  sumPzDau+=resdau->Pz();
1947  nKaons++;
1948  arrayDauLab[nFoundpKpi++]=indResDau;
1949  if(nFoundpKpi>3) return -1;
1950  }else if(TMath::Abs(pdgresdau)==211){
1951  sumPxDau+=resdau->Px();
1952  sumPyDau+=resdau->Py();
1953  sumPzDau+=resdau->Pz();
1954  nPions++;
1955  arrayDauLab[nFoundpKpi++]=indResDau;
1956  if(nFoundpKpi>3) return -1;
1957  }else if(TMath::Abs(pdgresdau)==2212){
1958  sumPxDau+=resdau->Px();
1959  sumPyDau+=resdau->Py();
1960  sumPzDau+=resdau->Pz();
1961  nProtons++;
1962  arrayDauLab[nFoundpKpi++]=indResDau;
1963  if(nFoundpKpi>3) return -1;
1964  }
1965  }
1966  }else{
1967  return -1;
1968  }
1969  }
1970  if(nPions!=1) return -1;
1971  if(nKaons!=1) return -1;
1972  if(nProtons!=1) return -1;
1973  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1974  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1975  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1976  if(nDau==3) return 1;
1977  else if(nDau==2){
1978  if(codeRes==313) return 2;
1979  else if(codeRes==2224) return 3;
1980  else if(codeRes==3124) return 4;
1981  }
1982  }
1983  return -1;
1984 
1985 }
1986 //____________________________________________________________________________
1987 Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1989 
1990  if(label<0) return -1;
1991  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1992  TParticle* part = mcEvent->Particle(label);
1993  if(!part || !mcPart) return -1;
1994  Int_t pdgD=part->GetPdgCode();
1995  if(TMath::Abs(pdgD)!=4122) return -1;
1996 
1997  Int_t nDau=mcPart->GetNDaughters();
1998  Int_t labelFirstDau = mcPart->GetDaughterFirst();
1999  Int_t nPions=0;
2000  Int_t nProtons=0;
2001  Double_t sumPxDau=0.;
2002  Double_t sumPyDau=0.;
2003  Double_t sumPzDau=0.;
2004  Int_t nFoundppi=0;
2005 
2006  Int_t codeV0=-1;
2007  if(nDau==2){
2008  for(Int_t iDau=0; iDau<nDau; iDau++){
2009  Int_t indDau = labelFirstDau+iDau;
2010  if(indDau<0) return -1;
2011  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
2012  TParticle* dau=mcEvent->Particle(indDau);
2013  if(!mcDau || !dau) return -1;
2014  Int_t pdgdau=dau->GetPdgCode();
2015  if(TMath::Abs(pdgdau)==211){
2016  nPions++;
2017  sumPxDau+=dau->Px();
2018  sumPyDau+=dau->Py();
2019  sumPzDau+=dau->Pz();
2020  arrayDauLab[nFoundppi++]=indDau;
2021  if(nFoundppi>3) return -1;
2022  }else if(TMath::Abs(pdgdau)==2212){
2023  nProtons++;
2024  sumPxDau+=dau->Px();
2025  sumPyDau+=dau->Py();
2026  sumPzDau+=dau->Pz();
2027  arrayDauLab[nFoundppi++]=indDau;
2028  if(nFoundppi>3) return -1;
2029  }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){
2030  codeV0=TMath::Abs(pdgdau);
2031  TParticle* v0=dau;
2032  AliMCParticle* mcV0=mcDau;
2033  if(codeV0==311){
2034  Int_t nK0Dau=mcDau->GetNDaughters();
2035  if(nK0Dau!=1) return -1;
2036  Int_t indK0s=mcDau->GetDaughterFirst();
2037  if(indK0s<0) return -1;
2038  v0=mcEvent->Particle(indK0s);
2039  mcV0=(AliMCParticle*)mcEvent->GetTrack(indK0s);
2040  if(!v0 || !mcV0) return -1;
2041  Int_t pdgK0sl=v0->GetPdgCode();
2042  codeV0=TMath::Abs(pdgK0sl);
2043  }
2044  Int_t nV0Dau=mcV0->GetNDaughters();
2045  if(nV0Dau!=2) return -1;
2046  Int_t indFirstV0Dau=mcV0->GetDaughterFirst();
2047  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
2048  Int_t indV0Dau=indFirstV0Dau+v0Dau;
2049  if(indV0Dau<0) return -1;
2050  TParticle* v0dau=mcEvent->Particle(indV0Dau);
2051  if(!v0dau) return -1;
2052  Int_t pdgv0dau=v0dau->GetPdgCode();
2053  if(TMath::Abs(pdgv0dau)==211){
2054  sumPxDau+=v0dau->Px();
2055  sumPyDau+=v0dau->Py();
2056  sumPzDau+=v0dau->Pz();
2057  nPions++;
2058  arrayDauLab[nFoundppi++]=indV0Dau;
2059  if(nFoundppi>3) return -1;
2060  }else if(TMath::Abs(pdgv0dau)==2212){
2061  sumPxDau+=v0dau->Px();
2062  sumPyDau+=v0dau->Py();
2063  sumPzDau+=v0dau->Pz();
2064  nProtons++;
2065  arrayDauLab[nFoundppi++]=indV0Dau;
2066  if(nFoundppi>3) return -1;
2067  }
2068  }
2069  }else{
2070  return -1;
2071  }
2072  }
2073  if(nPions!=2) return -1;
2074  if(nProtons!=1) return -1;
2075  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
2076  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
2077  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
2078  if(codeV0==310) return 1;
2079  else if(codeV0==3122) return 2;
2080  }
2081  return -1;
2082 
2083 }
2084 //____________________________________________________________________________
2085 Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
2087 
2088  Int_t pdgD=mcPart->GetPdgCode();
2089  if(TMath::Abs(pdgD)!=4122) return -1;
2090 
2091  Int_t nDau=mcPart->GetNDaughters();
2092  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
2093  Int_t nPions=0;
2094  Int_t nProtons=0;
2095  Double_t sumPxDau=0.;
2096  Double_t sumPyDau=0.;
2097  Double_t sumPzDau=0.;
2098  Int_t nFoundppi=0;
2099 
2100  Int_t codeV0=-1;
2101  if(nDau==2){
2102  for(Int_t iDau=0; iDau<nDau; iDau++){
2103  Int_t indDau = labelFirstDau+iDau;
2104  if(indDau<0) return -1;
2105  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
2106  if(!dau) return -1;
2107  Int_t pdgdau=dau->GetPdgCode();
2108  if(TMath::Abs(pdgdau)==211){
2109  nPions++;
2110  sumPxDau+=dau->Px();
2111  sumPyDau+=dau->Py();
2112  sumPzDau+=dau->Pz();
2113  arrayDauLab[nFoundppi++]=indDau;
2114  if(nFoundppi>3) return -1;
2115  }else if(TMath::Abs(pdgdau)==2212){
2116  nProtons++;
2117  sumPxDau+=dau->Px();
2118  sumPyDau+=dau->Py();
2119  sumPzDau+=dau->Pz();
2120  arrayDauLab[nFoundppi++]=indDau;
2121  if(nFoundppi>3) return -1;
2122  }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){
2123  codeV0=TMath::Abs(pdgdau);
2124  AliAODMCParticle* v0=dau;
2125  if(codeV0==311){
2126  Int_t nK0Dau=dau->GetNDaughters();
2127  if(nK0Dau!=1) return -1;
2128  Int_t indK0s=dau->GetDaughterLabel(0);
2129  if(indK0s<0) return -1;
2130  v0=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indK0s));
2131  if(!v0) return -1;
2132  Int_t pdgK0sl=v0->GetPdgCode();
2133  codeV0=TMath::Abs(pdgK0sl);
2134  }
2135  Int_t nV0Dau=v0->GetNDaughters();
2136  if(nV0Dau!=2) return -1;
2137  Int_t indFirstV0Dau=v0->GetDaughterLabel(0);
2138  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
2139  Int_t indV0Dau=indFirstV0Dau+v0Dau;
2140  if(indV0Dau<0) return -1;
2141  AliAODMCParticle* v0dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indV0Dau));
2142  if(!v0dau) return -1;
2143  Int_t pdgv0dau=v0dau->GetPdgCode();
2144  if(TMath::Abs(pdgv0dau)==211){
2145  sumPxDau+=v0dau->Px();
2146  sumPyDau+=v0dau->Py();
2147  sumPzDau+=v0dau->Pz();
2148  nPions++;
2149  arrayDauLab[nFoundppi++]=indV0Dau;
2150  if(nFoundppi>3) return -1;
2151  }else if(TMath::Abs(pdgv0dau)==2212){
2152  sumPxDau+=v0dau->Px();
2153  sumPyDau+=v0dau->Py();
2154  sumPzDau+=v0dau->Pz();
2155  nProtons++;
2156  arrayDauLab[nFoundppi++]=indV0Dau;
2157  if(nFoundppi>3) return -1;
2158  }
2159  }
2160  }else{
2161  return -1;
2162  }
2163  }
2164  if(nPions!=2) return -1;
2165  if(nProtons!=1) return -1;
2166  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2167  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2168  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2169  if(codeV0==310) return 1;
2170  else if(codeV0==3122) return 2;
2171  }
2172  return -1;
2173 
2174 }
2175 
2176 //__________________________________xic______________________________________
2177 Int_t AliVertexingHFUtils::CheckXicXipipiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2179 
2180  if(label<0) return -1;
2181  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
2182  TParticle* part = mcEvent->Particle(label);
2183  if(!part || !mcPart) return -1;
2184  Int_t pdgD=part->GetPdgCode();
2185  if(TMath::Abs(pdgD)!=4232) return -1;
2186 
2187  Int_t nDau=mcPart->GetNDaughters();
2188  if(nDau!=3) return -1;
2189 
2190  Int_t labelFirstDau = mcPart->GetDaughterFirst();
2191  Int_t nXi=0;
2192  Int_t nPions=0;
2193  Double_t sumPxDau=0.;
2194  Double_t sumPyDau=0.;
2195  Double_t sumPzDau=0.;
2196  Int_t nFoundXi=0;
2197 
2198  for(Int_t iDau=0; iDau<nDau; iDau++){
2199  Int_t indDau = labelFirstDau+iDau;
2200  if(indDau<0) return -1;
2201  TParticle* dau=mcEvent->Particle(indDau);
2202  if(!dau) return -1;
2203  Int_t pdgdau=dau->GetPdgCode();
2204  if(TMath::Abs(pdgdau)==3312){
2205  if(pdgD*pdgdau<0) return -1;
2206  sumPxDau+=dau->Px();
2207  sumPyDau+=dau->Py();
2208  sumPzDau+=dau->Pz();
2209  nXi++;
2210  arrayDauLab[nFoundXi++]=indDau;
2211 
2212  }
2213  if(TMath::Abs(pdgdau)==211){
2214  if(pdgD*pdgdau<0) return -1;
2215  nPions++;
2216  sumPxDau+=dau->Px();
2217  sumPyDau+=dau->Py();
2218  sumPzDau+=dau->Pz();
2219  arrayDauLab[nFoundXi++]=indDau;
2220  if(nFoundXi>3) return -1;
2221  }
2222  }
2223 
2224  if(nPions!=2) return -1;
2225  if(nXi!=1) return -1;
2226  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
2227  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
2228  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
2229  return 1;
2230 
2231 }
2232 
2233 //____________________________________________________________________________
2234 Int_t AliVertexingHFUtils::CheckBplusDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2236 
2237  if(label<0) return -1;
2238  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
2239  TParticle* part = mcEvent->Particle(label);
2240  if(!part || !mcPart) return -1;
2241  Int_t pdgD=part->GetPdgCode();
2242  if(TMath::Abs(pdgD)!=521) return -1;
2243 
2244  Int_t nDau=mcPart->GetNDaughters();
2245  if(nDau!=2) return -1;
2246 
2247  Int_t labelFirstDau = mcPart->GetDaughterFirst();
2248  Int_t nKaons=0;
2249  Int_t nPions=0;
2250  Double_t sumPxDau=0.;
2251  Double_t sumPyDau=0.;
2252  Double_t sumPzDau=0.;
2253  Int_t nFoundKpi=0;
2254 
2255  for(Int_t iDau=0; iDau<nDau; iDau++){
2256  Int_t indDau = labelFirstDau+iDau;
2257  if(indDau<0) return -1;
2258  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
2259  TParticle* dau=mcEvent->Particle(indDau);
2260  if(!mcDau || !dau) return -1;
2261  Int_t pdgdau=dau->GetPdgCode();
2262  if(TMath::Abs(pdgdau)==421){
2263  Int_t nResDau=mcDau->GetNDaughters();
2264  if(nResDau!=2) return -1;
2265  Int_t indFirstResDau=mcDau->GetDaughterFirst();
2266  for(Int_t resDau=0; resDau<2; resDau++){
2267  Int_t indResDau=indFirstResDau+resDau;
2268  if(indResDau<0) return -1;
2269  TParticle* resdau=mcEvent->Particle(indResDau);
2270  if(!resdau) return -1;
2271  Int_t pdgresdau=resdau->GetPdgCode();
2272  if(TMath::Abs(pdgresdau)==321){
2273  if(pdgD*pdgresdau<0) return -1;
2274  sumPxDau+=resdau->Px();
2275  sumPyDau+=resdau->Py();
2276  sumPzDau+=resdau->Pz();
2277  nKaons++;
2278  arrayDauLab[nFoundKpi++]=indResDau;
2279  if(nFoundKpi>3) return -1;
2280  }
2281  if(TMath::Abs(pdgresdau)==211){
2282  if(pdgD*pdgresdau>0) return -1;
2283  sumPxDau+=resdau->Px();
2284  sumPyDau+=resdau->Py();
2285  sumPzDau+=resdau->Pz();
2286  nPions++;
2287  arrayDauLab[nFoundKpi++]=indResDau;
2288  if(nFoundKpi>3) return -1;
2289  }
2290  }
2291  }else if(TMath::Abs(pdgdau)==211){
2292  if(pdgD*pdgdau<0) return -1;
2293  nPions++;
2294  sumPxDau+=dau->Px();
2295  sumPyDau+=dau->Py();
2296  sumPzDau+=dau->Pz();
2297  arrayDauLab[nFoundKpi++]=indDau;
2298  if(nFoundKpi>3) return -1;
2299  }
2300  }
2301 
2302  if(nPions!=2) return -1;
2303  if(nKaons!=1) return -1;
2304  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
2305  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
2306  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
2307  return 1;
2308 
2309 }
2310 //____________________________________________________________________________
2311 Int_t AliVertexingHFUtils::CheckBplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
2313 
2314  Int_t pdgD=mcPart->GetPdgCode();
2315  if(TMath::Abs(pdgD)!=521) return -1;
2316 
2317  Int_t nDau=mcPart->GetNDaughters();
2318  if(nDau!=2) return -1;
2319 
2320  Int_t labelFirstDau = mcPart->GetDaughterLabel(0);
2321  Int_t nKaons=0;
2322  Int_t nPions=0;
2323  Double_t sumPxDau=0.;
2324  Double_t sumPyDau=0.;
2325  Double_t sumPzDau=0.;
2326  Int_t nFoundKpi=0;
2327 
2328  for(Int_t iDau=0; iDau<nDau; iDau++){
2329  Int_t indDau = labelFirstDau+iDau;
2330  if(indDau<0) return -1;
2331  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
2332  if(!dau) return -1;
2333  Int_t pdgdau=dau->GetPdgCode();
2334  if(TMath::Abs(pdgdau)==421){
2335  Int_t nResDau=dau->GetNDaughters();
2336  if(nResDau!=2) return -1;
2337  Int_t indFirstResDau=dau->GetDaughterLabel(0);
2338  for(Int_t resDau=0; resDau<2; resDau++){
2339  Int_t indResDau=indFirstResDau+resDau;
2340  if(indResDau<0) return -1;
2341  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
2342  if(!resdau) return -1;
2343  Int_t pdgresdau=resdau->GetPdgCode();
2344  if(TMath::Abs(pdgresdau)==321){
2345  if(pdgD*pdgresdau<0) return -1;
2346  sumPxDau+=resdau->Px();
2347  sumPyDau+=resdau->Py();
2348  sumPzDau+=resdau->Pz();
2349  nKaons++;
2350  arrayDauLab[nFoundKpi++]=indResDau;
2351  if(nFoundKpi>3) return -1;
2352  }
2353  if(TMath::Abs(pdgresdau)==211){
2354  if(pdgD*pdgresdau>0) return -1;
2355  sumPxDau+=resdau->Px();
2356  sumPyDau+=resdau->Py();
2357  sumPzDau+=resdau->Pz();
2358  nPions++;
2359  arrayDauLab[nFoundKpi++]=indResDau;
2360  if(nFoundKpi>3) return -1;
2361  }
2362  }
2363  }else if(TMath::Abs(pdgdau)==211){
2364  if(pdgD*pdgdau<0) return -1;
2365  nPions++;
2366  sumPxDau+=dau->Px();
2367  sumPyDau+=dau->Py();
2368  sumPzDau+=dau->Pz();
2369  arrayDauLab[nFoundKpi++]=indDau;
2370  if(nFoundKpi>3) return -1;
2371  }
2372  }
2373 
2374  if(nPions!=2) return -1;
2375  if(nKaons!=1) return -1;
2376  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2377  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2378  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2379  return 1;
2380 
2381 }
2382 //________________________________________________________________________
2384  Double_t etaMin, Double_t etaMax,
2386  Int_t filtbit1, Int_t filtbit2,
2387  Int_t minMult){
2389 
2390  Int_t nTracks=aod->GetNumberOfTracks();
2391  Int_t nSelTracks=0;
2392 
2393  Double_t sumpt=0.;
2394  Double_t s00=0.;
2395  Double_t s01=0.;
2396  Double_t s11=0.;
2397  if(ptMin<0.) ptMin=0.;
2398 
2399  for(Int_t it=0; it<nTracks; it++) {
2400  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
2401  if(!tr) continue;
2402  Float_t eta = tr->Eta();
2403  Float_t pt = tr->Pt();
2404  Float_t phi = tr->Phi();
2405  if(eta<etaMin || eta>etaMax) continue;
2406  if(pt<ptMin || pt>ptMax) continue;
2407  Bool_t fb1 = tr->TestFilterBit(filtbit1);
2408  Bool_t fb2 = tr->TestFilterBit(filtbit2);
2409  Bool_t tpcRefit=tr->GetStatus() & AliAODTrack::kTPCrefit;
2410  if(filtbit1==1 && !tpcRefit) fb1=kFALSE;
2411  if(filtbit2==1 && !tpcRefit) fb2=kFALSE;
2412  if( !(fb1 || fb2) ) continue;
2413  Double_t px=pt*TMath::Cos(phi);
2414  Double_t py=pt*TMath::Sin(phi);
2415  s00 += (px * px)/pt;
2416  s01 += (py * px)/pt;
2417  s11 += (py * py)/pt;
2418  nSelTracks++;
2419  sumpt+=pt;
2420  }
2421 
2422  if(nSelTracks<minMult) return -0.5;
2423 
2424  if(sumpt>0.){
2425  s00/=sumpt;
2426  s01/=sumpt;
2427  s11/=sumpt;
2428  }else return -0.5;
2429 
2430  Double_t sphericity = -10;
2431  Double_t lambda1=((s00+s11)+TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
2432  Double_t lambda2=((s00+s11)-TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
2433  if(TMath::Abs(lambda2)<0.00001 && TMath::Abs(lambda1)<0.00001) sphericity=0;
2434  if(TMath::Abs(lambda1+lambda2)>0.000001) sphericity=2*TMath::Min(lambda1,lambda2)/(lambda1+lambda2);
2435  return sphericity;
2436 
2437 }
2438 
2439 //________________________________________________________________________
2441  Double_t &spherocity, Double_t &phiRef,
2442  Double_t etaMin, Double_t etaMax,
2444  Int_t filtbit1, Int_t filtbit2,
2445  Int_t minMult, Double_t phiStepSizeDeg,
2446  Int_t nTrksToSkip, Int_t* idToSkip
2447  ){
2449 
2450  Int_t nTracks=aod->GetNumberOfTracks();
2451  Int_t nSelTracks=0;
2452 
2453  Double_t* ptArr=new Double_t[nTracks];
2454  Double_t* phiArr=new Double_t[nTracks];
2455  Double_t sumpt=0.;
2456 
2457  for(Int_t it=0; it<nTracks; it++) {
2458  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
2459  if(!tr) continue;
2460  Float_t eta = tr->Eta();
2461  Float_t pt = tr->Pt();
2462  Float_t phi = tr->Phi();
2463  if(eta<etaMin || eta>etaMax) continue;
2464  if(pt<ptMin || pt>ptMax) continue;
2465  if(nTrksToSkip>0 && idToSkip){
2466  Int_t trid = (Int_t)tr->GetID();
2467  Bool_t keep=kTRUE;
2468  for(Int_t jt=0; jt<nTrksToSkip; jt++){
2469  if(trid==idToSkip[jt]) keep=kFALSE;
2470  }
2471  if(!keep) continue;
2472  }
2473  Bool_t fb1 = tr->TestFilterBit(filtbit1);
2474  Bool_t fb2 = tr->TestFilterBit(filtbit2);
2475  Bool_t tpcRefit=tr->GetStatus() & AliAODTrack::kTPCrefit;
2476  if(filtbit1==1 && !tpcRefit) fb1=kFALSE;
2477  if(filtbit2==1 && !tpcRefit) fb2=kFALSE;
2478  if( !(fb1 || fb2) ) continue;
2479  ptArr[nSelTracks]=pt;
2480  phiArr[nSelTracks]=phi;
2481  nSelTracks++;
2482  sumpt+=pt;
2483  }
2484 
2485  if(nSelTracks<minMult){spherocity = -0.5; return;}
2486 
2487  //Getting thrust
2488  spherocity=2.;
2489  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
2490  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
2491  Double_t nx=TMath::Cos(phistep);
2492  Double_t ny=TMath::Sin(phistep);
2493  Double_t numer=0.;
2494  for(Int_t j=0; j<nSelTracks; ++j){
2495  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
2496  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
2497  numer+=TMath::Abs(ny*pxA - nx*pyA);
2498  }
2499  Double_t pFull=numer*numer/(sumpt*sumpt);
2500  if(pFull<spherocity){
2501  spherocity=pFull; // minimization;
2502  phiRef=phistep;
2503  }
2504  }
2505 
2506  delete [] ptArr;
2507  delete [] phiArr;
2508 
2509  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
2510  return;
2511 
2512 }
2513 //________________________________________________________________________
2515  Double_t &spherocity, Double_t &phiRef,
2516  Double_t etaMin, Double_t etaMax,
2518  Int_t minMult, Double_t phiStepSizeDeg){
2519 
2521 
2522  Int_t nParticles=arrayMC->GetEntriesFast();
2523  Int_t nSelParticles=0;
2524 
2525  Double_t* ptArr=new Double_t[nParticles];
2526  Double_t* phiArr=new Double_t[nParticles];
2527  Double_t sumpt=0.;
2528 
2529  for(Int_t ip=0; ip<nParticles; ip++) {
2530  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(ip);
2531  if(!part) continue;
2532  Float_t eta = part->Eta();
2533  Float_t pt = part->Pt();
2534  Float_t phi = part->Phi();
2535  Int_t charge = part->Charge();
2536  Bool_t isPhysPrim = part->IsPhysicalPrimary();
2537  if(!isPhysPrim) continue;
2538  if(charge==0) continue;
2539  if(eta<etaMin || eta>etaMax) continue;
2540  if(pt<ptMin || pt>ptMax) continue;
2541 
2542  ptArr[nSelParticles]=pt;
2543  phiArr[nSelParticles]=phi;
2544  nSelParticles++;
2545  sumpt+=pt;
2546  }
2547 
2548  if(nSelParticles<minMult){spherocity = -0.5; return;}
2549 
2550  //Getting thrust
2551  spherocity=2.;
2552  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
2553  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
2554  Double_t nx=TMath::Cos(phistep);
2555  Double_t ny=TMath::Sin(phistep);
2556  Double_t numer=0.;
2557  for(Int_t j=0; j<nSelParticles; ++j){
2558  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
2559  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
2560  numer+=TMath::Abs(ny*pxA - nx*pyA);
2561  }
2562  Double_t pFull=numer*numer/(sumpt*sumpt);
2563  if(pFull<spherocity){
2564  spherocity=pFull; // minimization;
2565  phiRef=phistep;
2566  }
2567  }
2568 
2569  delete [] ptArr;
2570  delete [] phiArr;
2571 
2572  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
2573  return;
2574 
2575 }
2576 
2577 //________________________________________________________________________
2584 
2585  Int_t nBinOrig=hOrig->GetNbinsX();
2586  Int_t firstBinOrig=1;
2587  Int_t lastBinOrig=nBinOrig;
2588  Int_t nBinOrigUsed=nBinOrig;
2589  Int_t nBinFinal=nBinOrig/reb;
2590  if(firstUse>=1){
2591  firstBinOrig=firstUse;
2592  nBinFinal=(nBinOrig-firstUse+1)/reb;
2593  nBinOrigUsed=nBinFinal*reb;
2594  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
2595  }else{
2596  Int_t exc=nBinOrigUsed%reb;
2597  if(exc!=0){
2598  nBinOrigUsed-=exc;
2599  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
2600  }
2601  }
2602 
2603  printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
2604  Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
2605  Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
2606  TH1D* hRebin=new TH1D(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
2607  Int_t lastSummed=firstBinOrig-1;
2608  for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
2609  Float_t sum=0.;
2610  Float_t sume2=0.;
2611  for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
2612  sum+=hOrig->GetBinContent(lastSummed+1);
2613  sume2+=(hOrig->GetBinError(lastSummed+1)*hOrig->GetBinError(lastSummed+1));
2614  lastSummed++;
2615  }
2616  hRebin->SetBinContent(iBin,sum);
2617  hRebin->SetBinError(iBin,TMath::Sqrt(sume2));
2618  }
2619  return hRebin;
2620 }
2621 //________________________________________________________________________
2624 
2625  Int_t binmin=TMath::Max(1,hData->FindBin(hMC->GetXaxis()->GetXmin()));
2626  Bool_t found=kFALSE;
2627  Int_t binminD=-1;
2628  Int_t binminMC=-1;
2629  for(Int_t j=binmin; j<hData->GetNbinsX(); j++){
2630  if(found) break;
2631  for(Int_t k=1; k<hMC->GetNbinsX(); k++){
2632  Double_t delta=TMath::Abs(hMC->GetBinLowEdge(k)-hData->GetBinLowEdge(j));
2633  if(delta<0.0001){
2634  found=kTRUE;
2635  binminMC=k;
2636  binminD=j;
2637  }
2638  if(found) break;
2639  }
2640  }
2641  Int_t binmax=TMath::Min(hData->GetNbinsX(),hData->FindBin(hMC->GetXaxis()->GetXmax()*0.99999));
2642  found=kFALSE;
2643  Int_t binmaxD=-1;
2644  Int_t binmaxMC=-1;
2645  for(Int_t j=binmax; j>1; j--){
2646  if(found) break;
2647  for(Int_t k=hMC->GetNbinsX(); k>400; k--){
2648  Double_t delta=TMath::Abs(hMC->GetBinLowEdge(k+1)-hData->GetBinLowEdge(j+1));
2649  if(delta<0.0001){
2650  found=kTRUE;
2651  binmaxMC=k;
2652  binmaxD=j;
2653  }
2654  if(found) break;
2655  }
2656  }
2657 
2658  Double_t min=hData->GetBinLowEdge(binminD);
2659  Double_t max=hData->GetBinLowEdge(binmaxD)+hData->GetBinWidth(binmaxD);
2660  Double_t minMC=hMC->GetBinLowEdge(binminMC);
2661  Double_t maxMC=hMC->GetBinLowEdge(binmaxMC)+hMC->GetBinWidth(binmaxMC);
2662  Double_t width=hData->GetBinWidth(binminD);
2663  Double_t widthMC=hMC->GetBinWidth(binminMC);
2664 
2665  if(TMath::Abs(minMC-min)>0.0001*min || TMath::Abs(maxMC-max)>0.0001*max){
2666  printf("Cannot adapt range and rebin histo:\n");
2667  printf("Range for data histo: %f-%f GeV/c2 bins %d-%d width=%f\n",min,max,binminD,binmaxD,width);
2668  printf("Range for reflection histo: %f-%f GeV/c2 bins %d-%d width=%f\n",minMC,maxMC,binminMC,binmaxMC,widthMC);
2669  return 0x0;
2670  }
2671 
2672  Double_t rebin=width/widthMC;
2673  if(TMath::Abs(rebin-TMath::Nint(rebin))>0.001){
2674  printf("Cannot adapt histo: rebin %f issue, width MC = %f, width hData=%f (check=%f)\n",rebin,widthMC,width,TMath::Abs(rebin-TMath::Nint(rebin)));
2675  return 0x0;
2676  }
2677 
2678  Int_t nBinsNew=binmaxD-binminD+1;
2679  TH1 *hOut;
2680  TString stype=hMC->ClassName();
2681  if(stype.Contains("TH1F")){
2682  hOut=new TH1F(Form("%s-rebinned",hMC->GetName()),hMC->GetTitle(),nBinsNew,min,max);
2683  }else if(stype.Contains("TH1D")){
2684  hOut=new TH1D(Form("%s-rebinned",hMC->GetName()),hMC->GetTitle(),nBinsNew,min,max);
2685  }else{
2686  printf("Wrong type %s\n",stype.Data());
2687  return 0x0;
2688  }
2689 
2690  for(Int_t j=1; j<=hMC->GetNbinsX(); j++){
2691  Double_t m=hMC->GetBinCenter(j);
2692  Int_t binFin=hOut->FindBin(m);
2693  if(binFin>=1 && binFin<=nBinsNew){
2694  hOut->AddBinContent(binFin,hMC->GetBinContent(j));
2695  }
2696  }
2697  return hOut;
2698 }
2699 
2700 //___________________________________________________________________________________//
2701 //method that performs simultaneus fit of in-plane and out-of-plane inv-mass spectra
2702 ROOT::Fit::FitResult AliVertexingHFUtils::DoInPlaneOutOfPlaneSimultaneusFit(AliHFInvMassFitter *&massfitterInPlane, AliHFInvMassFitter *&massfitterOutOfPlane, TH1F* hMassInPlane, TH1F* hMassOutOfPlane, Double_t MinMass, Double_t MaxMass, Double_t massD, vector<UInt_t> commonpars) {
2703 
2704  cout << "\nIn-plane - out-of-plane simultaneus fit" << endl;
2705  cout << "\nIndependent prefits" << endl;
2706 
2707  //prefits to initialise parameters
2708  massfitterInPlane->SetUseLikelihoodFit();
2709  massfitterInPlane->SetInitialGaussianMean(massD);
2710  massfitterInPlane->SetInitialGaussianSigma(0.010);
2711  massfitterInPlane->MassFitter(kFALSE);
2712 
2713  massfitterOutOfPlane->SetUseLikelihoodFit();
2714  massfitterOutOfPlane->SetInitialGaussianMean(massD);
2715  massfitterOutOfPlane->SetInitialGaussianSigma(0.010);
2716  massfitterOutOfPlane->MassFitter(kFALSE);
2717 
2718  ROOT::Math::WrappedMultiTF1 wfInPlane(*(massfitterInPlane->GetMassFunc()),1);
2719  ROOT::Math::WrappedMultiTF1 wfOutOfPlane(*(massfitterOutOfPlane->GetMassFunc()),1);
2720 
2721  // set data options and ranges
2722  ROOT::Fit::DataOptions opt;
2723  ROOT::Fit::DataRange rangeMass; //same range for two functions
2724 
2725  rangeMass.SetRange(MinMass,MaxMass);
2726  ROOT::Fit::BinData dataInPlane(opt,rangeMass);
2727  ROOT::Fit::FillData(dataInPlane, hMassInPlane);
2728  ROOT::Fit::BinData dataOutOfPlane(opt,rangeMass);
2729  ROOT::Fit::FillData(dataOutOfPlane, hMassOutOfPlane);
2730 
2731  //define the 2 chi squares
2732  ROOT::Fit::Chi2Function chi2InPlane(dataInPlane, wfInPlane);
2733  ROOT::Fit::Chi2Function chi2OutOfPlane(dataOutOfPlane, wfOutOfPlane);
2734 
2735  //define the global chi square and get initial parameters from prefits
2736  const Int_t npars = massfitterInPlane->GetMassFunc()->GetNpar();
2737  const UInt_t ncommonpars = commonpars.size();
2738  Int_t nparsBkg = massfitterInPlane->GetBackgroundRecalcFunc()->GetNpar();
2739  Int_t nparsSgn = massfitterInPlane->GetSignalFunc()->GetNpar();
2740  Int_t nparsSecPeak = 0, nparsRefl = 0;
2741  if(massfitterInPlane->GetSecondPeakFunc())
2742  nparsSecPeak = massfitterInPlane->GetSecondPeakFunc()->GetNpar();
2743  if(massfitterInPlane->GetReflFunc())
2744  nparsRefl = massfitterInPlane->GetReflFunc()->GetNpar();
2745 
2746  GlobalInOutOfPlaneChi2 globalChi2(chi2InPlane, chi2OutOfPlane, npars, commonpars);
2747 
2748  vector<UInt_t>::iterator iter;
2749  vector<Double_t> initpars;
2750  for(Int_t iPar=0; iPar<2*npars; iPar++) {
2751  if(iPar<npars) { //in-plane
2752  initpars.push_back(massfitterInPlane->GetMassFunc()->GetParameter(iPar));
2753  }
2754  else { //out-of-plane
2755  iter = find(commonpars.begin(),commonpars.end(),iPar-npars);
2756  if(iter!=commonpars.end()) continue;
2757  else {
2758  initpars.push_back(massfitterOutOfPlane->GetMassFunc()->GetParameter(iPar-npars));
2759  }
2760  }
2761  }
2762 
2763  //define fitter and fit
2764  ROOT::Fit::Fitter simulfitter;
2765  simulfitter.Config().SetParamsSettings(npars*2-ncommonpars,initpars.data()); //set initial parameters from prefits
2766  if(nparsSecPeak>0) { //fix S/R
2767  simulfitter.Config().ParSettings(nparsBkg+nparsSgn+nparsSecPeak).Fix();
2768 
2769  iter = find(commonpars.begin(),commonpars.end(),npars+nparsBkg+nparsSgn+nparsSecPeak);
2770  if(iter==commonpars.end()) //if not included in common pars, need to be fixed also for out-of-plane func
2771  simulfitter.Config().ParSettings(npars+nparsBkg+nparsSgn+nparsSecPeak-ncommonpars).Fix();
2772  }
2773 
2774  simulfitter.Config().MinimizerOptions().SetPrintLevel(0);
2775  simulfitter.Config().SetMinimizer("Minuit2","Migrad");
2776  simulfitter.FitFCN(npars*2-ncommonpars,globalChi2,0,dataInPlane.Size()+dataOutOfPlane.Size(),kFALSE);
2777  ROOT::Fit::FitResult result = simulfitter.Result();
2778  cout << "\nSimultaneus fit" << endl;
2779  result.Print(cout);
2780 
2781  //Set new parameters to functions of mass fitters
2782  Int_t ncommonparsused = 0;
2783  for(Int_t iPar=0; iPar<npars; iPar++) {
2784  massfitterInPlane->GetMassFunc()->SetParameter(iPar,result.Parameter(iPar));
2785  massfitterInPlane->GetMassFunc()->SetParError(iPar,result.ParError(iPar));
2786  iter = find(commonpars.begin(),commonpars.end(),iPar);
2787  if(iter!=commonpars.end()) { //is common parameter
2788  massfitterOutOfPlane->GetMassFunc()->SetParameter(iPar,result.Parameter(iPar));
2789  massfitterOutOfPlane->GetMassFunc()->SetParError(iPar,result.ParError(iPar));
2790  ncommonparsused++;
2791  }
2792  else {
2793  massfitterOutOfPlane->GetMassFunc()->SetParameter(iPar,result.Parameter(iPar+npars-ncommonparsused));
2794  massfitterOutOfPlane->GetMassFunc()->SetParError(iPar,result.ParError(iPar+npars-ncommonparsused));
2795  }
2796 
2797  if(iPar < nparsBkg) { //bkg
2798  massfitterInPlane->GetBackgroundRecalcFunc()->SetParameter(iPar,massfitterInPlane->GetMassFunc()->GetParameter(iPar));
2799  massfitterInPlane->GetBackgroundRecalcFunc()->SetParError(iPar,massfitterInPlane->GetMassFunc()->GetParError(iPar));
2800  massfitterOutOfPlane->GetBackgroundRecalcFunc()->SetParameter(iPar,massfitterOutOfPlane->GetMassFunc()->GetParameter(iPar));
2801  massfitterOutOfPlane->GetBackgroundRecalcFunc()->SetParError(iPar,massfitterOutOfPlane->GetMassFunc()->GetParError(iPar));
2802  }
2803  else if(iPar >= nparsBkg && iPar < nparsBkg+nparsSgn){ //signal
2804  massfitterInPlane->GetSignalFunc()->SetParameter(iPar-nparsBkg,massfitterInPlane->GetMassFunc()->GetParameter(iPar));
2805  massfitterInPlane->GetSignalFunc()->SetParError(iPar-nparsBkg,massfitterInPlane->GetMassFunc()->GetParError(iPar));
2806  massfitterOutOfPlane->GetSignalFunc()->SetParameter(iPar-nparsBkg,massfitterOutOfPlane->GetMassFunc()->GetParameter(iPar));
2807  massfitterOutOfPlane->GetSignalFunc()->SetParError(iPar-nparsBkg,massfitterOutOfPlane->GetMassFunc()->GetParError(iPar));
2808  }
2809  else if(iPar >= nparsBkg && iPar < nparsBkg+nparsSgn){ //signal
2810  massfitterInPlane->GetSignalFunc()->SetParameter(iPar-nparsBkg,massfitterInPlane->GetMassFunc()->GetParameter(iPar));
2811  massfitterInPlane->GetSignalFunc()->SetParError(iPar-nparsBkg,massfitterInPlane->GetMassFunc()->GetParError(iPar));
2812  massfitterOutOfPlane->GetSignalFunc()->SetParameter(iPar-nparsBkg,massfitterOutOfPlane->GetMassFunc()->GetParameter(iPar));
2813  massfitterOutOfPlane->GetSignalFunc()->SetParError(iPar-nparsBkg,massfitterOutOfPlane->GetMassFunc()->GetParError(iPar));
2814  }
2815  }
2816 
2817  cout << "\n" << endl;
2818  return result;
2819 }
Int_t charge
static void GetSpherocity(AliAODEvent *aod, Double_t &spherocity, Double_t &phiRef, Double_t etaMin=-0.8, Double_t etaMax=0.8, Double_t ptMin=0.15, Double_t ptMax=10., Int_t filtbit1=256, Int_t filtbit2=512, Int_t minMult=3, Double_t phiStepSizeDeg=0.1, Int_t nTrksToSkip=0, Int_t *idToSkip=0x0)
Functions for event shape variables.
static Int_t CheckD0Decay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static void 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=2.5, Float_t sigmaRangeForBkg=4.5, Float_t minMass=0., Float_t maxMass=3., Int_t rebin=1)
Functions for computing average pt.
static TH1D * RebinHisto(TH1 *hOrig, Int_t reb, Int_t firstUse=-1)
Rebinning of invariant mass histograms.
double Double_t
Definition: External.C:58
Definition: External.C:236
static Int_t CheckDplusDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Int_t CheckDsDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Int_t CheckLcpKpiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
AliAODv0 * Getv0() const
Double_t ImpParXY() const
Int_t MassFitter(Bool_t draw=kTRUE)
static TString GetGenerator(Int_t label, AliAODMCHeader *header)
AliHFInvMassFitter class for the fit of invariant mass distribution of charm hadrons.
TCanvas * c
Definition: TestFitELoss.C:172
Double_t Pol(Double_t x) const
static Double_t GetFullEvResolLowLim(const TH1F *hSubEvCorr, Int_t k=1)
static Int_t CheckBplusDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
Double_t ptMin
static TH1 * AdaptTemplateRangeAndBinning(const TH1 *hRef, TH1 *hData, Double_t minFit, Double_t maxFit)
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Functions to check the decay tree.
TRandom * gRandom
void GetTrackPrimaryGenerator(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen)
Double_t fMinEtaForTracklets
sub-event resolution = sqrt(<cos[n(phiA-phiB)] >)
Double_t fSubRes
ratio of measured harmonic to event plane harmonic
static Int_t CheckDplusK0spiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Int_t CheckXicXipipiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Int_t CheckLcV0bachelorDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
Double_t fMaxEtaForTracklets
min eta for counting tracklets
static Double_t GetVZEROCEqualizedMultiplicity(AliAODEvent *ev)
static Double_t GetFullEvResolHighLim(const TH1F *hSubEvCorr, Int_t k=1)
static void GetGeneratedSpherocity(TClonesArray *arrayMC, Double_t &spherocity, Double_t &phiRef, Double_t etaMin=-0.8, Double_t etaMax=0.8, Double_t ptMin=0.15, Double_t ptMax=10., Int_t minMult=3, Double_t phiStepSizeDeg=0.1)
Double_t massD
static Double_t ResolK1(Double_t x)
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
int Int_t
Definition: External.C:63
Double_t GetFullEvResol() const
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
const Double_t ptmax
static Int_t CheckDplusKKpiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Bool_t CheckT0TriggerFired(AliAODEvent *aodEv)
Functions for processing trigger information.
Definition: External.C:212
AliAODTrack * GetBachelor() const
static Double_t GetSubEvResolLowLim(const TH1F *hSubEvCorr)
const Double_t ptmin
Bool_t HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header, TClonesArray *arrayMC)
Bool_t IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header, TClonesArray *arrayMC)
static Int_t CheckDstarDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
short Short_t
Definition: External.C:23
static Double_t GetBeautyMotherPt(TClonesArray *arrayMC, AliAODMCParticle *mcPart)
static void ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance, Double_t &errsignificance)
Significance calculator.
static Double_t GetCorrectedNtracklets(TProfile *estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult)
static ROOT::Fit::FitResult DoInPlaneOutOfPlaneSimultaneusFit(AliHFInvMassFitter *&massfitterInPlane, AliHFInvMassFitter *&massfitterOutOfPlane, TH1F *hMassInPlane, TH1F *hMassOutOfPlane, Double_t MinMass, Double_t MaxMass, Double_t massD, vector< UInt_t > commonpars)
Double_t GetSubEvResol() const
Bool_t IsTrackInjected(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC)
void SetInitialGaussianSigma(Double_t sigma)
void SetInitialGaussianMean(Double_t mean)
Double_t minMass
Double_t FindChi() const
Int_t rebin
static Int_t GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray *arrayMC, Double_t mineta, Double_t maxeta)
static Double_t GetVZEROAEqualizedMultiplicity(AliAODEvent *ev)
Utilities for V0 multiplicity checks.
static Double_t GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partDp)
Functions for computing true impact parameter of D meson.
Double_t maxMass
bool Bool_t
Definition: External.C:53
static Int_t CheckDsK0sKDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
Double_t ptMax
static Double_t GetSubEvResolHighLim(const TH1F *hSubEvCorr)
static Double_t GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partDp)
static Double_t GetSphericity(AliAODEvent *aod, Double_t etaMin=-0.8, Double_t etaMax=0.8, Double_t ptMin=0.15, Double_t ptMax=10., Int_t filtbit1=256, Int_t filtbit2=512, Int_t minMult=3)
Definition: External.C:196
static Int_t GetGeneratedMultiplicityInEtaRange(TClonesArray *arrayMC, Double_t mineta, Double_t maxeta)
static Int_t GetGeneratedPrimariesInEtaRange(TClonesArray *arrayMC, Double_t mineta, Double_t maxeta)
Class with functions useful for different D2H analyses //.