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