AliPhysics  95775ff (95775ff)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliHFsubtractBFDcuts.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-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 
25 #include "AliHFsubtractBFDcuts.h"
26 
27 #include <vector>
28 
29 #include "TNamed.h"
30 #include "THnSparse.h"
31 #include "TH1F.h"
32 #include "TH3F.h"
33 #include "TClonesArray.h"
34 #include "TString.h"
35 #include "TObjString.h"
36 #include "TMath.h"
37 
39 #include "AliAODMCHeader.h"
40 #include "AliAODMCParticle.h"
41 #include "AliAODEvent.h"
42 #include "AliAODVertex.h"
43 #include "AliESDVertex.h"
44 #include "AliVTrack.h"
45 #include "AliAODTrack.h"
46 #include "AliVertexerTracks.h"
47 #include "AliExternalTrackParam.h"
48 #include "AliNeutralTrackParam.h"
49 
53 
55  : TNamed()
56  , fIsMC(kTRUE)
57  , fCheckAcceptance(kTRUE)
58  , fResolveResonances(kTRUE)
59  , fTHnGenStep(0x0)
60  , fTHnData(0x0)
61  , fTHnMC(0x0)
62  , fMCarray(0x0)
63  , fAODtracks(0x0)
64  , fLabCand(-1)
65  , fPtCand(0.)
66  , fNprongs((UInt_t)-1)
67  , fNprongsInAcc((UInt_t)-1)
68  , fFoundElectron(kFALSE)
69  , fMotherPt(-1.)
70  , fGenerateDecayList(kFALSE) // DON'T ACTIVATE, DOESN'T MERGE
71  , fDecayProngs()
72  , fDecayStrList(0x0)
73  , fQAhists(0x0)
74 {
76 }
77 
78 AliHFsubtractBFDcuts::AliHFsubtractBFDcuts(const char* name, const char* title)
79  : TNamed(name,title)
80  , fIsMC(kTRUE)
81  , fCheckAcceptance(kTRUE)
82  , fResolveResonances(kTRUE)
83  , fTHnGenStep(0x0)
84  , fTHnData(0x0)
85  , fTHnMC(0x0)
86  , fMCarray(0x0)
87  , fAODtracks(0x0)
88  , fLabCand(-1)
89  , fPtCand(0.)
90  , fNprongs((UInt_t)-1)
91  , fNprongsInAcc((UInt_t)-1)
92  , fFoundElectron(kFALSE)
93  , fMotherPt(-1.)
94  , fGenerateDecayList(kFALSE) // DON'T ACTIVATE, DOESN'T MERGE
95  , fDecayProngs()
96  , fDecayStrList(0x0)
97  , fQAhists(0x0)
98 {
100  fDecayStrList = new TList();
101  fDecayStrList->SetOwner();
102  fQAhists = new TList();
103  fQAhists->SetOwner();
104 }
105 
107  : TNamed(c.GetName(),c.GetTitle())
108  , fIsMC(c.fIsMC)
109  , fCheckAcceptance(c.fCheckAcceptance)
110  , fResolveResonances(c.fResolveResonances)
111  , fTHnGenStep(c.fTHnGenStep)
112  , fTHnData(c.fTHnData)
113  , fTHnMC(c.fTHnMC)
114  , fMCarray(c.fMCarray)
115  , fAODtracks(c.fAODtracks)
116  , fLabCand(c.fLabCand)
117  , fPtCand(c.fPtCand)
118  , fNprongs(c.fNprongs)
119  , fNprongsInAcc(c.fNprongsInAcc)
120  , fFoundElectron(c.fFoundElectron)
121  , fMotherPt(c.fMotherPt)
122  , fGenerateDecayList(c.fGenerateDecayList)
123  , fDecayProngs(c.fDecayProngs)
124  , fDecayStrList(c.fDecayStrList) // FIXME: TList copy contructor not implemented
125  , fQAhists(c.fQAhists)
126 {
128 }
129 
131 {
133  return c;
134 }
135 
137  AliDebug(3, "TODO: implement me!");
138 }
139 
141  // mass, pt, normLXY, cosPointXY, normL, cosPoint, LXY, L
142  Int_t dimAxes[8]={500 ,96 ,200 ,100 , 200,100 , 100 ,100 };
143  Double_t min[8]={ 1.7, 0., 0., 0.99, 0., 0.99, 0. , 0. };
144  Double_t max[8]={ 2.2,24.,100., 1. ,100., 1. , 1.0, 1.0};
145  fTHnData=new THnSparseF("fCutsDataFD","fCutsDataFD",8,dimAxes,min,max);
146  fTHnData->GetAxis(0)->SetName("mass");
147  fTHnData->GetAxis(0)->SetTitle("Mass (K,#pi) (GeV/#it{c^{2}})");
148  fTHnData->GetAxis(1)->SetName("pt");
149  fTHnData->GetAxis(1)->SetTitle("#it{p}_{T} (GeV/#it{c})");
150  fTHnData->GetAxis(2)->SetName("NormDecLengthXY");
151  fTHnData->GetAxis(2)->SetTitle("Normalized XY decay length");
152  fTHnData->GetAxis(3)->SetName("CosPointXY");
153  fTHnData->GetAxis(3)->SetTitle("Cos#theta_{point}^{XY}");
154  fTHnData->GetAxis(4)->SetName("NormDecLength");
155  fTHnData->GetAxis(4)->SetTitle("Normalized decay length");
156  fTHnData->GetAxis(5)->SetName("CosPoint");
157  fTHnData->GetAxis(5)->SetTitle("Cos#theta_{point}");
158  fTHnData->GetAxis(6)->SetName("DecLengthXY");
159  fTHnData->GetAxis(6)->SetTitle("XY decay length (cm)");
160  fTHnData->GetAxis(7)->SetName("DecLength");
161  fTHnData->GetAxis(7)->SetTitle("Decay length (cm)");
162 
163  // pt, normLXY, cosPointXY, #prongs, mother pt, normL, cosPoint, LXY, L
164  Int_t dimAxesMC[10]={96 ,200 ,100 ,20 ,24 ,200 ,100 ,100 ,100 ,2 };
165  Double_t minMC[10]={ 0., 0., 0.99, 0., 0., 0., 0.99, 0. , 0. ,0.};
166  Double_t maxMC[10]={24.,100., 1. ,20.,24.,100., 1. , 1.0, 1.0,2.};
167  fTHnMC=new THnSparseF("fCutsMCFD","fCutsMCFD",10,dimAxesMC,minMC,maxMC);
168  fTHnMC->GetAxis(0)->SetName("pt");
169  fTHnMC->GetAxis(0)->SetTitle("#it{p}_{T} (GeV/#it{c})");
170  fTHnMC->GetAxis(1)->SetName("NormDecLengthXY");
171  fTHnMC->GetAxis(1)->SetTitle("Normalized XY decay length");
172  fTHnMC->GetAxis(2)->SetName("CosPointXY");
173  fTHnMC->GetAxis(2)->SetTitle("Cos#theta_{point}^{XY}");
174  fTHnMC->GetAxis(3)->SetName("nProngs");
175  fTHnMC->GetAxis(3)->SetTitle("Number of Decay Prongs");
176  fTHnMC->GetAxis(4)->SetName("Mother_pt");
177  fTHnMC->GetAxis(4)->SetTitle("Mother #it{p}_{T} (GeV/#it{c})");
178  fTHnMC->GetAxis(5)->SetName("NormDecLength");
179  fTHnMC->GetAxis(5)->SetTitle("Normalized decay length");
180  fTHnMC->GetAxis(6)->SetName("CosPoint");
181  fTHnMC->GetAxis(6)->SetTitle("Cos#theta_{point}");
182  fTHnMC->GetAxis(7)->SetName("DecLengthXY");
183  fTHnMC->GetAxis(7)->SetTitle("XY decay length (cm)");
184  fTHnMC->GetAxis(8)->SetName("DecLength");
185  fTHnMC->GetAxis(8)->SetTitle("Decay length (cm)");
186  fTHnMC->GetAxis(9)->SetName("ContainsElectron");
187  fTHnMC->GetAxis(9)->SetTitle("ContainsElectron");
188 
189  Int_t dimAxesGen[6]={96 ,20 ,24 ,100 ,100 ,2 };
190  Double_t minGen[6]={ 0., 0., 0., 0., 0.,0.};
191  Double_t maxGen[6]={24.,20.,24., 1., 1.,2.};
192  fTHnGenStep=new THnSparseF("fPtMCGenStep","fPtMCGenStep",6,dimAxesGen,minGen,maxGen);
193  fTHnGenStep->GetAxis(0)->SetName("pt");
194  fTHnGenStep->GetAxis(0)->SetTitle("#it{p}_{T} (GeV/#it{c})");
195  fTHnGenStep->GetAxis(1)->SetName("nProngs");
196  fTHnGenStep->GetAxis(1)->SetTitle("Number of Decay Prongs");
197  fTHnGenStep->GetAxis(2)->SetName("Mother_pt");
198  fTHnGenStep->GetAxis(2)->SetTitle("Mother #it{p}_{T} (GeV/#it{c})");
199  fTHnGenStep->GetAxis(3)->SetName("DecLengthXY");
200  fTHnGenStep->GetAxis(3)->SetTitle("XY decay length (cm)");
201  fTHnGenStep->GetAxis(4)->SetName("DecLength");
202  fTHnGenStep->GetAxis(4)->SetTitle("Decay length (cm)");
203  fTHnGenStep->GetAxis(5)->SetName("ContainsElectron");
204  fTHnGenStep->GetAxis(5)->SetTitle("ContainsElectron");
205 
206  fQAhists->Add(new TH1F("hRapidityDist" , "All Particles;y;Counts (a.u.)" ,800,-4., 4. )); // 0
207  fQAhists->Add(new TH1F("hRapidityDistStable" , "All Stable Particles;y;Counts (a.u.)" ,800,-4., 4. )); // 1
208  fQAhists->Add(new TH2F("hDecayLengthB" , ";D meson #it{p}_{T};Decay Length B Meson (cm);counts (a.u.)" , 24, 0.,24. , 200,0.,2.0 )); // 2
209  fQAhists->Add(new TH2F("hDecayLengthBxy" , ";D meson #it{p}_{T};Decay Length B Meson in XY direction (cm);counts (a.u.)" , 24, 0.,24. , 200,0.,2.0 )); // 3
210  fQAhists->Add(new TH3F("hDCABhypothesis" , "All Particles;D meson #it{p}_{T};Distance of Closest Approach (cm);Impact parameter of D0 x track (cm^{2})" , 24, 0.,24. , 100,0.,1.0 ,100,0.,0.001)); // 4
211  fQAhists->Add(new TH3F("hDCABprongs" , "Only Real B Decay Prongs;D meson #it{p}_{T};Distance of Closest Approach (cm);Impact parameter of D0 x track (cm^{2})" , 24, 0.,24. , 100,0.,1.0 ,100,0.,0.001)); // 5
212  fQAhists->Add(new TH3F("hVtxDistBprongs" , "Only Real B Decay Prongs;D meson #it{p}_{T};Distance to Primary Vertex (cm);Error Distance to Primary Vertex (cm)" , 24, 0.,24. , 100,0.,1.0 ,100,0.,1.0 )); // 6
213  fQAhists->Add(new TH3F("hVtxDistXYBprongs" , "Only Real B Decay Prongs;D meson #it{p}_{T};Distance to Primary Vertex in XY (cm);Error Distance to Primary Vertex in XY (cm)", 24, 0.,24. , 100,0.,1.0 ,100,0.,1.0 )); // 7
214  fQAhists->Add(new TH3F("hVtxDistBhypothesis" , "All Particles;D meson #it{p}_{T};Distance to Primary Vertex (cm);Error Distance to Primary Vertex (cm)" , 24, 0.,24. , 100,0.,1.0 ,100,0.,1.0 )); // 8
215  fQAhists->Add(new TH3F("hVtxDistXYBhypthesis" , "All Particles;D meson #it{p}_{T};Distance to Primary Vertex in XY (cm);Error Distance to Primary Vertex in XY (cm)" , 24, 0.,24. , 100,0.,1.0 ,100,0.,1.0 )); // 9
216  fQAhists->Add(new TH2F("hNormVtxDistBprongs" , "Only Real B Decay Prongs;D meson #it{p}_{T};Distance to Primary Vertex (cm);counts (a.u.)" , 24, 0.,24. , 100,0.,1.0 )); // 10
217  fQAhists->Add(new TH2F("hNormVtxDistXYBprongs" , "Only Real B Decay Prongs;D meson #it{p}_{T};Distance to Primary Vertex in XY (cm);counts (a.u.)" , 24, 0.,24. , 100,0.,1.0 )); // 11
218  fQAhists->Add(new TH2F("hNormVtxDistBhypothesis" , "All Particles;D meson #it{p}_{T};Distance to Primary Vertex (cm);counts (a.u.)" , 24, 0.,24. , 100,0.,1.0 )); // 12
219  fQAhists->Add(new TH2F("hNormVtxDistXYBhypthesis", "All Particles;D meson #it{p}_{T};Distance to Primary Vertex in XY (cm);counts (a.u.)" , 24, 0.,24. , 100,0.,1.0 )); // 13
220  fQAhists->Add(new TH2F("hVtxPrecision" , ";#it{d}_{z}^{Vertex} (cm);#it{d}_{xyz}^{Vertex} (cm);counts (a.u.)" ,100, 0., 0.05,1000,0.,0.05 )); // 14
221  fQAhists->Add(new TH2F("hDecayLengthPrecision" , ";D meson #it{p}_{T};Decay Length Residual (cm);counts (a.u.)" , 24, 0.,24. ,1000,0.,0.05 )); // 15
222  fQAhists->Add(new TH2F("hDecayLengthXYPrecision" , ";D meson #it{p}_{T};XY decay Length Residual(cm);counts (a.u.)" , 24, 0.,24. ,1000,0.,0.05 )); // 16
223  fQAhists->ls();
224 
225  return;
226 }
227 
228 void AliHFsubtractBFDcuts::FillGenStep(AliAODMCParticle* dzeropart,Double_t pt/*=-1.*/,Double_t weight/*=1.*/,TClonesArray* mcArray/*=0x0*/,AliAODMCHeader* mcHeader/*=0x0*/){
229  fMCarray=mcArray;
230  if (pt<0) {
231  pt=dzeropart->Pt();
232  }
233  if (fIsMC && fMCarray) {
234  fNprongs=0;
235  fNprongsInAcc=0;
236  fFoundElectron=kFALSE;
237  fDecayChain=kFALSE; // TODO: use this value
238  fMotherPt=pt;
239  fLabCand=dzeropart->GetLabel();
240  fPtCand=dzeropart->Pt();
241  Double_t vtxDist[] = { mcHeader->GetVtxX()-((AliAODMCParticle*)mcArray->At(dzeropart->GetFirstDaughter()))->Xv(),
242  mcHeader->GetVtxY()-((AliAODMCParticle*)mcArray->At(dzeropart->GetFirstDaughter()))->Yv(),
243  mcHeader->GetVtxZ()-((AliAODMCParticle*)mcArray->At(dzeropart->GetFirstDaughter()))->Zv() };
244  Double_t decayLength = TMath::Sqrt(vtxDist[0]*vtxDist[0]+vtxDist[1]*vtxDist[1]+vtxDist[2]*vtxDist[2]);
245  Double_t decayLengthXY = TMath::Sqrt(vtxDist[0]*vtxDist[0]+vtxDist[1]*vtxDist[1]);
246  if (!AnalyseDecay(fGenerateDecayList, kTRUE)) {
247  AliDebug(3, "Error during the decay type determination!");
248  }
249  Double_t entry[] = {pt,(Double_t)fNprongsInAcc,fMotherPt,decayLength,decayLengthXY,(Double_t)fFoundElectron};
250  fTHnGenStep->Fill(entry,weight);
251  // y distribution of all particles
252  for (Int_t i=0; i<fMCarray->GetEntriesFast(); ++i) {
253  Double_t y=((AliAODMCParticle*)fMCarray->UncheckedAt(i))->Y();
254  ((TH1F*)fQAhists->At(0))->Fill(y); // all particles
255  if (IsStable(i)) ((TH1F*)fQAhists->At(1))->Fill(y); // all stable particles
256  }
257  }
258  else {
259  Double_t entry[] = {pt,0.,-1., 0., 0.};
260  fTHnGenStep->Fill(entry,weight);
261  }
262  fMCarray=0x0;
263  return;
264 }
265 
266 void AliHFsubtractBFDcuts::FillSparses(AliAODRecoDecayHF2Prong* dzerocand,Int_t isSelected,Double_t pt,Double_t massD0,Double_t massD0bar,Double_t weight,TClonesArray* mcArray,AliAODEvent* aodEvent, AliAODMCHeader* mcHeader){
267  fMCarray=mcArray;
268  if (aodEvent) {
269  fAODtracks=aodEvent->GetTracks();
270  fBkG =aodEvent->GetMagneticField();
271  fPriVtx =aodEvent->GetPrimaryVertex();
272  }
273  if(isSelected<=0||isSelected>3){
274  Printf("isSelected = %d", isSelected);
275  return;
276  }
277  if(massD0<0){
278  dzerocand->InvMassD0(massD0,massD0bar);
279  Printf("mass D0 =%f, massD0bar = %f ", massD0, massD0bar);
280  }
281  if(pt<0){
282  pt=dzerocand->Pt();
283  }
284  Double_t normalDecayLengXY=dzerocand->NormalizedDecayLengthXY();
285  Double_t cptangXY=dzerocand->CosPointingAngleXY();
286  Double_t normalDecayLeng=dzerocand->NormalizedDecayLength();
287  Double_t cptang=dzerocand->CosPointingAngle();
288  Double_t decayLengXY=dzerocand->DecayLengthXY();
289  Double_t decayLeng=dzerocand->DecayLength();
290  // mass, pt, normLXY, cosPointXY, normL, cosPoint, LXY, L
291  Double_t pointData[]={massD0,pt,normalDecayLengXY,cptangXY,normalDecayLeng,cptang,decayLengXY,decayLeng};
292 
293  if(isSelected==1||isSelected==3){
294  fTHnData->Fill(pointData,weight);
295  }
296  if(isSelected==2||isSelected==3){
297  pointData[0]=massD0bar;
298  fTHnData->Fill(pointData,weight);
299  }
300 
301  if(fIsMC && fMCarray) {
302  fNprongs=0;
303  fNprongsInAcc=0;
304  fDecayChain=kFALSE; // TODO: use this value
305  fFoundElectron=kFALSE;
306  fMotherPt=pt;
307  fD0Cand=dzerocand;
308  if (fD0Cand) fD0CandParam=new AliNeutralTrackParam(fD0Cand);
309  if (!GetCandidateLabel() || !AnalyseDecay(kFALSE, kFALSE)) {
310  AliDebug(3, "Error during the decay type determination!");
311  }
312  // pt, normLXY, cosPointXY, #prongs, mother pt, normL, cosPoint, LXY, L
313  Double_t pointMC[]={pt,normalDecayLengXY,cptangXY,(Double_t)fNprongsInAcc,fMotherPt,normalDecayLeng,cptang,decayLengXY,decayLeng,(Double_t)fFoundElectron};
314  fTHnMC->Fill(pointMC, weight);
315 
316  if (mcHeader) {
317  Double_t dist[] = { fPriVtx->GetX()-mcHeader->GetVtxX(),
318  fPriVtx->GetY()-mcHeader->GetVtxY(),
319  fPriVtx->GetZ()-mcHeader->GetVtxZ() };
320  Int_t labDau0=((AliAODTrack*)fD0Cand->GetDaughter(0))->GetLabel();
321  if (labDau0<0) {
322  labDau0 *= -1.;
323  AliDebug(1, "Negative Daughter label");
324  }
325  Double_t decLengMC3D[] = { mcHeader->GetVtxX()-((AliAODMCParticle*)mcArray->At(labDau0))->Xv(),
326  mcHeader->GetVtxY()-((AliAODMCParticle*)mcArray->At(labDau0))->Yv(),
327  mcHeader->GetVtxZ()-((AliAODMCParticle*)mcArray->At(labDau0))->Zv() };
328  Double_t decayLengMC = TMath::Sqrt(decLengMC3D[0]*decLengMC3D[0]+decLengMC3D[1]*decLengMC3D[1]+decLengMC3D[2]*decLengMC3D[2]);
329  Double_t decayLengXYMC = TMath::Sqrt(decLengMC3D[0]*decLengMC3D[0]+decLengMC3D[1]*decLengMC3D[1]);
330 
331 
332  ((TH2F*)fQAhists->At(14))->Fill(TMath::Sqrt(dist[2]*dist[2]),
333  TMath::Sqrt(dist[0]*dist[0]+dist[1]*dist[1]+dist[2]*dist[2]));
334  ((TH2F*)fQAhists->At(15))->Fill(pt, decayLengMC-decayLeng);
335  ((TH2F*)fQAhists->At(16))->Fill(pt, decayLengXYMC-decayLengXY);
336  }
337  }
338  fMCarray=0x0;
339  fAODtracks=0x0;
340  fBkG=0.;
341  fD0Cand=0x0;
342  return;
343 }
344 
346  if (fD0Cand) {
347  Int_t labDau0=((AliAODTrack*)fD0Cand->GetDaughter(0))->GetLabel();
348  if (labDau0<0) {
349  labDau0 *= -1.;
350  AliDebug(1, "Negative Daughter label");
351  }
352  AliAODMCParticle* firstDau=(AliAODMCParticle*)fMCarray->UncheckedAt(labDau0);
353  fLabCand = firstDau->GetMother();
354  return kTRUE;
355  }
356  else {
357  AliDebug(3, "Could not obtain the label of the candidate");
358  return kFALSE;
359  }
360  return kFALSE;
361 }
362 
364  AliAODMCParticle* cand=(AliAODMCParticle*)fMCarray->UncheckedAt(fLabCand);
365  fLabMother = cand->GetMother();
366  AliAODMCParticle* mother = (AliAODMCParticle*)fMCarray->UncheckedAt(fLabMother);
367  Int_t pdgMother = mother->GetPdgCode();
368  if (pdgMother<0) pdgMother*=-1; // treat particles and anti-particles the same way
369  if (pdgMother==4 || pdgMother==2212 || pdgMother==2112) { // prompt production
370  fNprongs=1;
371  fNprongsInAcc=1;
372  return kTRUE;
373  }
374  if ((pdgMother%1000)/100==4 || (pdgMother%10000)/1000==4) {
375  // chained decay of charmed hadrons, using recursion to resolve it
376  fDecayChain=kTRUE;
378  return AnalyseDecay(generateString, mcOnly);
379  }
380  if ((pdgMother%1000)/100!=5 && (pdgMother%10000)/1000!=5) {
381  AliDebug(3, "Found strange decay, expected the mother to be a beauty hadron!");
382  fNprongs=0;
383  fNprongsInAcc=0;
384  return kFALSE;
385  }
386  CountProngs(fLabMother, fLabCand, generateString, mcOnly); // count the prongs
387 
388  if (generateString) {
389  // Store the decay
390  std::sort(fDecayProngs.begin(), fDecayProngs.end());
391  TString decayStr = "";
392  for (ULong64_t i=0; i<fDecayProngs.size(); ++i) {
393  decayStr = (i==0) ? Form("%d", fDecayProngs[i]) : Form("%s_%d", decayStr.Data(), fDecayProngs[i]);
394  }
395  //decayStr = Form("%s__%d_%d_%d", decayStr.Data(), fNprongs, fNprongsInAcc, fDecayProngs.size());
396  TObjString* str = new TObjString(decayStr);
397  if (!fDecayStrList->FindObject(str)) fDecayStrList->Add(str); // only allow unique entries
398  fDecayProngs.clear();
399  }
400 
401  if (mcOnly) {
402  Double_t decayLengthB; // Decay length B meson
403  Double_t decayLengthBxy; // Decay length B meson (xy-plane)
404  Double_t originB[3] = {0.,0.,0.};
405  if(!mother->XvYvZv(originB)) AliDebug(3, "Couldn't determine MC origin of the beauty hadron");
406  Double_t originD[3] = {0.,0.,0.};
407  if(!cand->XvYvZv(originD)) AliDebug(3, "Couldn't determine MC origin of the charmed hadron");
408  decayLengthBxy = TMath::Sqrt((originB[0]-originD[0])*(originB[0]-originD[0])+
409  (originB[1]-originD[1])*(originB[1]-originD[1]));
410  decayLengthB = TMath::Sqrt(decayLengthBxy*decayLengthBxy+(originB[2]-originD[2])*(originB[2]-originD[2]));
411  ((TH2F*)fQAhists->At(2))->Fill(fPtCand,decayLengthB);
412  ((TH2F*)fQAhists->At(3))->Fill(fPtCand,decayLengthBxy);
413  }
414  else {
415  for (Int_t iAODtrack=0; iAODtrack<fAODtracks->GetEntriesFast(); ++iAODtrack) {
416  CheckBhypothesis(iAODtrack, kFALSE);
417  }
418  }
419 
420  fMotherPt=mother->Pt();
421  return kTRUE;
422 }
423 
424 void AliHFsubtractBFDcuts::CountProngs(Int_t labCurrMother, Int_t labCurrExcl,
425  Bool_t generateString, Bool_t mcOnly) {
426  for (Int_t iMCParticle=0; iMCParticle<fMCarray->GetEntriesFast(); ++iMCParticle) {
427  if (iMCParticle!=labCurrExcl) {
428  if (((AliAODMCParticle*)fMCarray->UncheckedAt(iMCParticle))->GetMother()==labCurrMother) {
429  if (!fResolveResonances || IsStable(iMCParticle)) {
430  if (generateString) fDecayProngs.push_back(((AliAODMCParticle*)fMCarray->UncheckedAt(iMCParticle))->GetPdgCode());
431  ++fNprongs;
432  if (!mcOnly) {
433  for (Int_t iAODtrack=0; iAODtrack<fAODtracks->GetEntriesFast(); ++iAODtrack) {
434  AliAODTrack* aodTrack=(AliAODTrack*)fAODtracks->UncheckedAt(iAODtrack);
435  if (aodTrack->GetLabel()==iMCParticle) CheckBhypothesis(iAODtrack, kTRUE);
436  }
437  }
438  if (!fCheckAcceptance || IsInAcceptance(iMCParticle)) {
439  ++fNprongsInAcc;
440  }
441  }
442  else CountProngs(iMCParticle, -1, generateString, mcOnly);
443  }
444  }
445  else {
446  ++fNprongs; // candidate is only counted as a single prong
447  ++fNprongsInAcc;
448  if (generateString) fDecayProngs.push_back(((AliAODMCParticle*)fMCarray->UncheckedAt(iMCParticle))->GetPdgCode());
449  }
450  }
451 }
452 
454  const Int_t stablePartPdgs[] = { 11, 13, 211, 321, 2212, 12, 14, 22, 111, 130 };
455  const Int_t nStablePartPdgs = sizeof(stablePartPdgs)/sizeof(Int_t);
456  AliAODMCParticle* prong = (AliAODMCParticle*)fMCarray->UncheckedAt(labProng);
457  Int_t pdgProng = prong->GetPdgCode();
458  if (pdgProng<0) pdgProng*=-1; // treat particles and anti-particles the same way
459  for (Int_t iPdg=0; iPdg<nStablePartPdgs; ++iPdg) {
460  if (pdgProng == 11) fFoundElectron = kTRUE;
461  if (stablePartPdgs[iPdg] == pdgProng) return kTRUE;
462  }
463  return kFALSE;
464 }
465 
467  AliDebug(1, "AliHFsubtractBFDcuts::IsInAcceptance(...) hasn't been implemented yet, prong");
468  Double_t eta = ((AliAODMCParticle*)fMCarray->UncheckedAt(labProng))->Eta();
469  Double_t pt = ((AliAODMCParticle*)fMCarray->UncheckedAt(labProng))->Pt();
470  Short_t charge = ((AliAODMCParticle*)fMCarray->UncheckedAt(labProng))->Charge();
471  if ((pt>0.15) && (eta>-0.9) && (eta<0.9) && (charge!=0)) return kTRUE;
472  return kFALSE;
473 }
474 
476  AliExternalTrackParam* t = new AliExternalTrackParam();
477  t->CopyFromVTrack((AliVTrack*)fAODtracks->UncheckedAt(iAODtrack));
478 
479  TObjArray* tracks = new TObjArray(2);
480  tracks->AddAt(t, 0);
481  tracks->AddAt(fD0CandParam,1);
482 
483  AliAODVertex* bVtx = RecBvtx(tracks);
484  if(!bVtx) {
485  AliDebug(3, "Couldn't reconstruct B meson vertex!");
486  delete t; t=0x0;
487  return kFALSE;
488  }
489 
490  Double_t vtxDist = (fPriVtx) ? bVtx->DistanceToVertex(fPriVtx) : -1.0 ;
491  Double_t vtxDistXY = (fPriVtx) ? bVtx->DistanceXYToVertex(fPriVtx) : -2.e10;
492  Double_t vtxDistErr = (fPriVtx) ? bVtx->ErrorDistanceToVertex(fPriVtx) : -1.0 ;
493  Double_t vtxDistErrXY = (fPriVtx) ? bVtx->ErrorDistanceXYToVertex(fPriVtx) : -2.e10;
494 
495  const Double_t maxD = 1.;
496 
497  // Propagate candidates to secondary vertex
498  Double_t dz[2],cov[3];
499  t->PropagateToDCA(bVtx,fBkG,maxD,dz,cov);
500  fD0CandParam->PropagateToDCA(bVtx,fBkG,maxD,dz,cov);
501 
502  // Impact parameters
503  t->PropagateToDCA(fPriVtx,fBkG,maxD,dz,cov);
504  Double_t d0Track=dz[0];
505  Double_t d0TrackErr=TMath::Sqrt(cov[0]);
506  fD0CandParam->PropagateToDCA(fPriVtx,fBkG,maxD,dz,cov);
507  Double_t d0D0Cand=dz[0];
508  Double_t d0D0CandErr=TMath::Sqrt(cov[0]);
509 
510  // distance of closest approach of the D0 and the track
511  Double_t xDCAtrack, xDCAD0;
512  Double_t dcaB=t->GetDCA(fD0CandParam,fBkG,xDCAtrack,xDCAD0);
513  if (!Bprong) {
514  ((TH3F*)fQAhists->At(4))->Fill(fPtCand,dcaB,d0D0Cand*d0Track);
515  ((TH3F*)fQAhists->At(6))->Fill(fPtCand,vtxDist,vtxDistErr);
516  ((TH3F*)fQAhists->At(7))->Fill(fPtCand,vtxDistXY,vtxDistErrXY);
517  ((TH2F*)fQAhists->At(10))->Fill(fPtCand,vtxDist/vtxDistErr);
518  ((TH2F*)fQAhists->At(11))->Fill(fPtCand,vtxDistXY/vtxDistErrXY);
519  }
520  else {
521  ((TH3F*)fQAhists->At(5))->Fill(fPtCand,dcaB,d0D0Cand*d0Track);
522  ((TH3F*)fQAhists->At(8))->Fill(fPtCand,vtxDist,vtxDistErr);
523  ((TH3F*)fQAhists->At(9))->Fill(fPtCand,vtxDistXY,vtxDistErrXY);
524  ((TH2F*)fQAhists->At(12))->Fill(fPtCand,vtxDist/vtxDistErr);
525  ((TH2F*)fQAhists->At(13))->Fill(fPtCand,vtxDistXY/vtxDistErrXY);
526  }
527 
528  delete tracks; tracks=0x0;
529  delete bVtx ; bVtx =0x0;
530  delete t ; t =0x0;
531  return kTRUE;
532 }
533 
534 AliAODVertex* AliHFsubtractBFDcuts::RecBvtx(TObjArray* tracks) const {
535  AliESDVertex* vtxESD=0x0;
536  AliVertexerTracks* vertexer = new AliVertexerTracks(fBkG);
537  vertexer->SetVtxStart((AliESDVertex*)fPriVtx);
538  vtxESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(tracks);
539  delete vertexer; vertexer=0x0;
540  if(!vtxESD || vtxESD->GetNContributors()!=tracks->GetEntriesFast()) {
541  AliDebug(3, "Couldn't reconstruct B Meson vertex");
542  delete vtxESD; vtxESD=0x0;
543  return 0x0;
544  }
545  Double_t rVtxSq=vtxESD->GetX()*vtxESD->GetX()+vtxESD->GetY()*vtxESD->GetY();
546  if(rVtxSq>8.){
547  // vertex outside beam pipe, reject candidate to avoid propagation through material
548  delete vtxESD; vtxESD=0x0;
549  return 0x0;
550  }
551  // convert to AliAODVertex
552  Double_t pos[3],cov[6],chi2perNDF;
553  vtxESD->GetXYZ(pos);
554  vtxESD->GetCovMatrix(cov);
555  chi2perNDF=vtxESD->GetChi2toNDF();
556  delete vtxESD; vtxESD=0x0;
557  return new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,0);
558 }
Int_t charge
Double_t NormalizedDecayLengthXY() const
Double_t NormalizedDecayLength() const
return jsonbuilder str().c_str()
Bool_t AnalyseDecay(Bool_t generateString, Bool_t mcOnly)
double Double_t
Definition: External.C:58
Definition: External.C:260
Bool_t fIsMC
Method to check Whether the current D0 candidate and the track originate from a B decay...
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
THnSparseF * fTHnMC
! THnSparse for cut variables (MC at PID level, w/o mass axis)y
AliAODVertex * RecBvtx(TObjArray *tracks) const
Is that prong within the fiducial acceptance.
THnSparseF * fTHnGenStep
flag resolve resonances in during the prong determination
Double_t fPtCand
Label of the candidate D0 (charmed hadron in case of a chained decay)
Class for storing and handling D0 meson candidates properties // for estimating the feed-down fractio...
std::vector< Int_t > fDecayProngs
Generate the list containig strings with all PDG codes of the decay prongs.
TCanvas * c
Definition: TestFitELoss.C:172
Double_t CosPointingAngleXY() const
AliAODRecoDecayHF2Prong * fD0Cand
Magnetic field (z-direction) in units of kG.
AliNeutralTrackParam * fD0CandParam
Pointer to the D0 candidate from reconstruction.
Int_t fLabMother
pT of the candidate (from MC track, not following decay chain)
Bool_t fCheckAcceptance
flag for MC/Data
TClonesArray * fAODtracks
! TClonesArray holding the AliAODTracks of the event to be processed
AliAODVertex * fPriVtx
! Primary AOD vertex
int Int_t
Definition: External.C:63
TList * fQAhists
! List with QA histograms
unsigned int UInt_t
Definition: External.C:33
TList * fDecayStrList
PDG codes of the daughters separated.
Bool_t fFoundElectron
Number of prongs, counting only the particles within acceptance.
TClonesArray * fMCarray
Event specific variables.
UInt_t fNprongsInAcc
Number of prongs, counting the first charmed hadron as one particle (simulation cuts can lead to loss...
UInt_t fNprongs
Label of the mother of the candidate D0 (or charmed hadron)
void CountProngs(Int_t labCurrMother, Int_t labCurrExcl, Bool_t generateString, Bool_t mConly)
check in which decay process a particle was created
Double_t fMotherPt
Chained decay of charmed hadrons.
short Short_t
Definition: External.C:23
Double_t DecayLengthXY() const
Bool_t fResolveResonances
flag for checking whether the decay prongs are within acceptance
Bool_t IsStable(Int_t labProng)
counting the prongs of labCurrMother, labCurrExcl is assumed to be a stable particle ...
Bool_t fDecayChain
Does the B meson decay contain an electron?
Bool_t IsInAcceptance(Int_t labProng) const
Is that prong a stable particle?
void FillSparses(AliAODRecoDecayHF2Prong *dzeroPart, Int_t isSelected, Double_t pt=-1, Double_t massD0=-1, Double_t massD0bar=-1, Double_t weight=1., TClonesArray *mcArray=0x0, AliAODEvent *aodEvent=0x0, AliAODMCHeader *mcHeader=0x0)
AliHFsubtractBFDcuts operator=(const AliHFsubtractBFDcuts &c)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
void FillGenStep(AliAODMCParticle *dzeroMC, Double_t pt=-1, Double_t weight=1., TClonesArray *mcArray=0x0, AliAODMCHeader *mcHeader=0x0)
bool Bool_t
Definition: External.C:53
Int_t fLabCand
Pointer to an AliNeutralTrackParam of the D0 candidata for DCA calculation.
Double_t CosPointingAngle() const
Double_t DecayLength() const
THnSparseF * fTHnData
! THnSparse for cut variables (data, with inv mass axis), first axis is always mass ...
Bool_t CheckBhypothesis(Int_t iAODtrack, Bool_t Bprong)
Reconstruct a secondary vertex with the supplied tracks.