AliPhysics  63e47e1 (63e47e1)
AliAnalysisTaskSEDStarCharmFraction.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 
16 //*************************************************************************
17 // Class AliAnalysisTaskSEDStarCharmFraction
18 // AliAnalysisTask for the extraction of the fraction of prompt charm for D*
19 // using the impact parameter to the primary vertex
20 //
21 // Author: Jasper van der Maarel <J.vanderMaarel@uu.nl>
22 //*************************************************************************
23 
24 #include <TSystem.h>
25 #include <TParticle.h>
26 #include <TH1D.h>
27 #include <THnSparse.h>
28 #include <TTree.h>
29 #include "TROOT.h"
30 #include <TDatabasePDG.h>
31 #include <TParameter.h>
32 #include <AliAnalysisDataSlot.h>
33 #include <AliAnalysisDataContainer.h>
35 #include "AliMCEvent.h"
36 #include "AliAnalysisManager.h"
37 #include "AliAODMCHeader.h"
38 #include "AliAODHandler.h"
39 #include "AliLog.h"
40 #include "AliAODVertex.h"
41 #include "AliAODRecoDecay.h"
42 #include "AliAODRecoDecayHF.h"
43 #include "AliAODRecoCascadeHF.h"
45 #include "AliAnalysisVertexingHF.h"
46 #include "AliESDtrack.h"
47 #include "AliAODMCParticle.h"
49 #include "AliAODEvent.h"
50 #include "AliVertexerTracks.h"
51 #include "AliNeutralTrackParam.h"
53 
55 
58  fCuts(0),
59  fCounter(0),
60  fReadMC(kFALSE),
61  fSkipHijing(kTRUE),
62  fSingleSideband(kFALSE),
63  fImpParCut(0),
64  fPDGMDStarD0(0),
65  fNPtBins(0),
66  fNEvents(0),
67  fTreeCandidate(0),
68  fListCandidate(0),
69  fListSignal(0),
70  fListSignalPrompt(0),
71  fListSignalFromB(0),
72  fListBackground(0),
73  fIsSideband(kFALSE),
74  fIsPeak(kFALSE),
75  fIsSignal(kFALSE),
76  fIsSignalPrompt(kFALSE),
77  fIsSignalFromB(kFALSE),
78  fIsBackground(kFALSE),
79  fNewPrimVtx(0),
80  fDStarVtx(0),
81  fMagneticField(0),
82  fPtForTree(0),
83  fInvMassForTree(0),
84  fImpParForTree(0),
85  fTrueImpParForTree(0),
86  fTypeForTree(0),
87  fSignalTypeForTree(0)
88 { // Default constructor
89 }
90 
92  AliAnalysisTaskSE(name),
93  fCuts(0),
94  fCounter(0),
95  fReadMC(kFALSE),
96  fSkipHijing(kTRUE),
97  fSingleSideband(kFALSE),
98  fImpParCut(0),
99  fPDGMDStarD0(0),
100  fNPtBins(0),
101  fNEvents(0),
102  fTreeCandidate(0),
103  fListCandidate(0),
104  fListSignal(0),
105  fListSignalPrompt(0),
106  fListSignalFromB(0),
107  fListBackground(0),
108  fIsSideband(kFALSE),
109  fIsPeak(kFALSE),
110  fIsSignal(kFALSE),
111  fIsSignalPrompt(kFALSE),
112  fIsSignalFromB(kFALSE),
113  fIsBackground(kFALSE),
114  fNewPrimVtx(0),
115  fDStarVtx(0),
116  fMagneticField(0),
117  fPtForTree(0),
118  fInvMassForTree(0),
119  fImpParForTree(0),
120  fTrueImpParForTree(0),
121  fTypeForTree(0),
122  fSignalTypeForTree(0)
123 { // Constructor
124 
125  if (fCuts) {
126  delete fCuts;
127  fCuts = 0;
128  }
129 
130  fCuts = new AliRDHFCutsDStartoKpipi(*cuts);
132 
133  DefineOutput(1, TH1D::Class()); // Event counter
134  DefineOutput(2, TList::Class()); // Candidate
135  DefineOutput(3, TList::Class()); // Signal
136  DefineOutput(4, TList::Class()); // SignalPrompt
137  DefineOutput(5, TList::Class()); // SignalFromB
138  DefineOutput(6, TList::Class()); // Background
139  DefineOutput(7, AliRDHFCutsDStartoKpipi::Class()); // Cuts
140  DefineOutput(8, AliNormalizationCounter::Class()); // Normalization
141  DefineOutput(9, TTree::Class()); // Tree for candidates
142 
143  for (Int_t i=0;i<30;i++) {
144  fPeakCut[i] = 0; fSidebandCut[i] = 0; fSidebandWindow[i] = 0;
145  }
146 
147  //DefineOutput(5, AliNormalizationCounter::Class());
148 }
149 
151 { // Destructor
152  if (fCuts) {
153  delete fCuts;
154  fCuts = 0;
155  }
156 
157  if (fNewPrimVtx) {
158  delete fNewPrimVtx;
159  fNewPrimVtx = 0;
160  }
161 
162  if (fDStarVtx) {
163  delete fDStarVtx;
164  fDStarVtx = 0;
165  }
166 }
167 
169 { // Initialization
170 
172  const char* name = GetOutputSlot(7)->GetContainer()->GetName();
173  copyCuts->SetName(name);
174 
175  PostData(7, copyCuts);
176 
177  fPDGMDStarD0 = TDatabasePDG::Instance()->GetParticle(413)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass();
178 }
179 
181 { // Create histograms
182  fNEvents = new TH1D(GetOutputSlot(1)->GetContainer()->GetName(), "Counter", 21, -0.5, 20.5);
183  fNEvents->GetXaxis()->SetBinLabel(1, "Events");
184  fNEvents->GetXaxis()->SetBinLabel(2, "Magnetic field");
185  fNEvents->GetXaxis()->SetBinLabel(3, "Selected");
186  fNEvents->GetXaxis()->SetBinLabel(4, "Primary vertex");
187  fNEvents->GetXaxis()->SetBinLabel(5, "Candidates");
188  fNEvents->GetXaxis()->SetBinLabel(6, "Passed track cuts");
189  fNEvents->GetXaxis()->SetBinLabel(7, "Passed fiducial acc cuts");
190  fNEvents->GetXaxis()->SetBinLabel(8, "Calc D* vtx");
191  fNEvents->GetXaxis()->SetBinLabel(9, "Pass daught imppar cut");
192  fNEvents->GetXaxis()->SetBinLabel(10, "Real D* (MC)");
193  fNEvents->GetXaxis()->SetBinLabel(11, "Skip from Hijing (MC)");
194  fListCandidate = new TList();
195  fListCandidate->SetOwner();
196  fListCandidate->SetName("listCandidate");
198 
199  fListSignal = new TList();
200  fListSignal->SetOwner();
201  fListSignal->SetName("listSignal");
202  if (fReadMC) {
204  }
205 
206  fListSignalPrompt = new TList();
207  fListSignalPrompt->SetOwner();
208  fListSignalPrompt->SetName("listSignalPrompt");
209  if (fReadMC) {
211  }
212 
213  fListSignalFromB = new TList();
214  fListSignalFromB->SetOwner();
215  fListSignalFromB->SetName("listSignalFromB");
216  if (fReadMC) {
218  }
219 
220  fListBackground = new TList();
221  fListBackground->SetOwner();
222  fListBackground->SetName("listBackground");
223  if (fReadMC) {
225  }
226 
227  fCounter = new AliNormalizationCounter(GetOutputSlot(8)->GetContainer()->GetName());
228  fCounter->Init();
229 
230  fTreeCandidate = new TTree(GetOutputSlot(9)->GetContainer()->GetName(), GetOutputSlot(9)->GetContainer()->GetName());
231  fTreeCandidate->Branch("pt", &fPtForTree, "pt/D");
232  fTreeCandidate->Branch("M", &fInvMassForTree, "M/D");
233  fTreeCandidate->Branch("d0", &fImpParForTree, "d0/D");
234  fTreeCandidate->Branch("trued0", &fTrueImpParForTree, "trued0/D");
235  fTreeCandidate->Branch("type", &fTypeForTree, "type/S");
236  fTreeCandidate->Branch("stype", &fSignalTypeForTree, "stype/S");
237 
238  PostData(1, fNEvents);
239  PostData(2, fListCandidate);
240  PostData(3, fListSignal);
241  PostData(4, fListSignalPrompt);
242  PostData(5, fListSignalFromB);
243  PostData(6, fListBackground);
244  PostData(8, fCounter);
245  PostData(9, fTreeCandidate);
246 }
247 
249 { // Fill a TList with histograms
250 
251  TString listName = list->GetName();
252  listName.ReplaceAll("list", "");
253  listName.ToLower();
254  if (listName.EqualTo("signalfromb")) {
255  listName = "signal from B";
256  }
257  if (listName.EqualTo("signalprompt")) {
258  listName = "signal prompt";
259  }
260 
261  TString regions[3] = {"", "Peak", "Sideband"}; // Invariant mass regions: all, peak region, sideband region
262 
263  // DeltaInvMass
264  TH1D* hDeltaInvMass = new TH1D("DeltaInvMass", Form("D*-D^{0} invariant mass, %s; #DeltaM [GeV/c^{2}]; Entries", listName.Data()), 700, 0.13, 0.2);
265  hDeltaInvMass->Sumw2();
266  hDeltaInvMass->SetLineColor(4);
267  hDeltaInvMass->SetMarkerColor(4);
268  hDeltaInvMass->SetMarkerStyle(20);
269  hDeltaInvMass->SetMarkerSize(0.6);
270  list->Add(hDeltaInvMass);
271 
272 
273 
274  // InvMassD0
275  TH1D* hInvMassD0 = new TH1D("InvMassD0", Form("D^{0} invariant mass, %s; M [GeV/c^{2}]; Entries", listName.Data()), 600, 1.6, 2.2);
276  hInvMassD0->Sumw2();
277  hInvMassD0->SetLineColor(4);
278  hInvMassD0->SetMarkerColor(4);
279  hInvMassD0->SetMarkerStyle(20);
280  hInvMassD0->SetMarkerSize(0.6);
281  list->Add(hInvMassD0);
282 
283  for (Int_t r=0;r<3;r++) {
284  TString regionTitle(regions[r]);
285  if (!regionTitle.EqualTo("")) {
286  regionTitle.ToLower();
287  regionTitle.Prepend(", ");
288  regionTitle.Append(" region");
289  }
290 
291  // ImpPar
292  TH1D* hImpPar = new TH1D(Form("ImpPar%s", regions[r].Data()), Form("D* impact parameter%s, %s; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data()), 1000, -1000., 1000.);
293  hImpPar->Sumw2();
294  hImpPar->SetLineColor(4);
295  hImpPar->SetMarkerColor(4);
296  hImpPar->SetMarkerStyle(20);
297  hImpPar->SetMarkerSize(0.6);
298  list->Add(hImpPar);
299 
300  // ImpParSoftPi
301  TH1D* hImpParSoftPi = new TH1D(Form("ImpParSoftPi%s", regions[r].Data()), Form("#pi_{s} impact parameter%s, %s; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data()), 1000, -1000., 1000.);
302  hImpParSoftPi->Sumw2();
303  hImpParSoftPi->SetLineColor(4);
304  hImpParSoftPi->SetMarkerColor(4);
305  hImpParSoftPi->SetMarkerStyle(20);
306  hImpParSoftPi->SetMarkerSize(0.6);
307  list->Add(hImpParSoftPi);
308 
309  // ImpParD0
310  TH1D* hImpParD0 = new TH1D(Form("ImpParD0%s", regions[r].Data()), Form("D^{0} impact parameter%s, %s; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data()), 1000, -1000., 1000.);
311  hImpParD0->Sumw2();
312  hImpParD0->SetLineColor(4);
313  hImpParD0->SetMarkerColor(4);
314  hImpParD0->SetMarkerStyle(20);
315  hImpParD0->SetMarkerSize(0.6);
316  list->Add(hImpParD0);
317 
318  //ImpParD0K
319  TH1D* hImpParD0K = new TH1D(Form("ImpParD0K%s", regions[r].Data()), Form("K (from D^{0}) impact parameter%s, %s; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data()), 1000, -1000., 1000.);
320  hImpParD0K->Sumw2();
321  hImpParD0K->SetLineColor(4);
322  hImpParD0K->SetMarkerColor(4);
323  hImpParD0K->SetMarkerStyle(20);
324  hImpParD0K->SetMarkerSize(0.6);
325  list->Add(hImpParD0K);
326 
327  //ImpParD0K2
328  TH1D* hImpParD0K2 = new TH1D(Form("ImpParD0K2%s", regions[r].Data()), Form("K (from D^{0}) impact parameter (no recalc prim vtx)%s, %s; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data()), 1000, -1000., 1000.);
329  hImpParD0K2->Sumw2();
330  hImpParD0K2->SetLineColor(4);
331  hImpParD0K2->SetMarkerColor(4);
332  hImpParD0K2->SetMarkerStyle(20);
333  hImpParD0K2->SetMarkerSize(0.6);
334  list->Add(hImpParD0K2);
335 
336  //ImpParD0Pi
337  TH1D* hImpParD0Pi = new TH1D(Form("ImpParD0Pi%s", regions[r].Data()), Form("#pi (from D^{0}) impact parameter%s, %s; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data()), 1000, -1000., 1000.);
338  hImpParD0Pi->Sumw2();
339  hImpParD0Pi->SetLineColor(4);
340  hImpParD0Pi->SetMarkerColor(4);
341  hImpParD0Pi->SetMarkerStyle(20);
342  hImpParD0Pi->SetMarkerSize(0.6);
343  list->Add(hImpParD0Pi);
344 
345  //ImpParD0Pi2
346  TH1D* hImpParD0Pi2 = new TH1D(Form("ImpParD0Pi2%s", regions[r].Data()), Form("#pi (from D^{0}) impact parameter (no recalc prim vtx)%s, %s; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data()), 1000, -1000., 1000.);
347  hImpParD0Pi2->Sumw2();
348  hImpParD0Pi2->SetLineColor(4);
349  hImpParD0Pi2->SetMarkerColor(4);
350  hImpParD0Pi2->SetMarkerStyle(20);
351  hImpParD0Pi2->SetMarkerSize(0.6);
352  list->Add(hImpParD0Pi2);
353  }
354 
355  Float_t *ptLimits = fCuts->GetPtBinLimits();
356 
357  for (Int_t i=0;i<fNPtBins;i++) {
358  Float_t ptLow = ptLimits[i];
359  Float_t ptHigh = ptLimits[(i+1)];
360 
361  TH1D* hDeltaInvMassLoop = new TH1D(Form("DeltaInvMass_%d", i), Form("D*-D^{0} invariant mass, %s, %.1f < p_{T} < %.1f GeV/c; #DeltaM [GeV/c^{2}]; Entries", listName.Data(), ptLow, ptHigh), 700, 0.13, 0.2);
362  hDeltaInvMassLoop->Sumw2();
363  hDeltaInvMassLoop->SetLineColor(4);
364  hDeltaInvMassLoop->SetMarkerColor(4);
365  hDeltaInvMassLoop->SetMarkerStyle(20);
366  hDeltaInvMassLoop->SetMarkerSize(0.6);
367  list->Add(hDeltaInvMassLoop);
368 
369 
370 
371  TH1D* hInvMassD0Loop = new TH1D(Form("InvMassD0_%d", i), Form("D^{0} invariant mass, %s, %.1f < p_{T} < %.1f GeV/c; M [GeV/c^{2}]; Entries", listName.Data(), ptLow, ptHigh), 600, 1.6, 2.2);
372  hInvMassD0Loop->Sumw2();
373  hInvMassD0Loop->SetLineColor(4);
374  hInvMassD0Loop->SetMarkerColor(4);
375  hInvMassD0Loop->SetMarkerStyle(20);
376  hInvMassD0Loop->SetMarkerSize(0.6);
377  list->Add(hInvMassD0Loop);
378 
379  for (Int_t r=0;r<3;r++) {
380  TString regionTitle(regions[r]);
381  if (!regionTitle.EqualTo("")) {
382  regionTitle.ToLower();
383  regionTitle.Prepend(", ");
384  regionTitle.Append(" region");
385  }
386 
387  //ImpPar
388  TH1D* hImpPar = new TH1D(Form("ImpPar%s_%d", regions[r].Data(), i), Form("D* impact parameter%s, %s, %.1f < p_{T} < %.1f GeV/c; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data(), ptLow, ptHigh), 1000, -1000., 1000.);
389  hImpPar->Sumw2();
390  hImpPar->SetLineColor(4);
391  hImpPar->SetMarkerColor(4);
392  hImpPar->SetMarkerStyle(20);
393  hImpPar->SetMarkerSize(0.6);
394  list->Add(hImpPar);
395 
396  //ImpParSoftPi
397  TH1D* hImpParSoftPi = new TH1D(Form("ImpParSoftPi%s_%d", regions[r].Data(), i), Form("#pi_{s} impact parameter%s, %s, %.1f < p_{T} < %.1f GeV/c; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data(), ptLow, ptHigh), 1000, -1000., 1000.);
398  hImpParSoftPi->Sumw2();
399  hImpParSoftPi->SetLineColor(4);
400  hImpParSoftPi->SetMarkerColor(4);
401  hImpParSoftPi->SetMarkerStyle(20);
402  hImpParSoftPi->SetMarkerSize(0.6);
403  list->Add(hImpParSoftPi);
404 
405  //ImpParD0
406  TH1D* hImpParD0 = new TH1D(Form("ImpParD0%s_%d", regions[r].Data(), i), Form("D^{0} impact parameter%s, %s, %.1f < p_{T} < %.1f GeV/c; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data(), ptLow, ptHigh), 1000, -1000., 1000.);
407  hImpParD0->Sumw2();
408  hImpParD0->SetLineColor(4);
409  hImpParD0->SetMarkerColor(4);
410  hImpParD0->SetMarkerStyle(20);
411  hImpParD0->SetMarkerSize(0.6);
412  list->Add(hImpParD0);
413 
414  //ImpParD0K
415  TH1D* hImpParD0K = new TH1D(Form("ImpParD0K%s_%d", regions[r].Data(), i), Form("K (from D^{0}) impact parameter%s, %s, %.1f < p_{T} < %.1f GeV/c; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data(), ptLow, ptHigh), 1000, -1000., 1000.);
416  hImpParD0K->Sumw2();
417  hImpParD0K->SetLineColor(4);
418  hImpParD0K->SetMarkerColor(4);
419  hImpParD0K->SetMarkerStyle(20);
420  hImpParD0K->SetMarkerSize(0.6);
421  list->Add(hImpParD0K);
422 
423  //ImpParD0K2
424  TH1D* hImpParD0K2 = new TH1D(Form("ImpParD0K2%s_%d", regions[r].Data(), i), Form("K (from D^{0}) impact parameter (no recalc prim vtx)%s, %s, %.1f < p_{T} < %.1f GeV/c; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data(), ptLow, ptHigh), 1000, -1000., 1000.);
425  hImpParD0K2->Sumw2();
426  hImpParD0K2->SetLineColor(4);
427  hImpParD0K2->SetMarkerColor(4);
428  hImpParD0K2->SetMarkerStyle(20);
429  hImpParD0K2->SetMarkerSize(0.6);
430  list->Add(hImpParD0K2);
431 
432  //ImpParD0Pi
433  TH1D* hImpParD0Pi = new TH1D(Form("ImpParD0Pi%s_%d", regions[r].Data(), i), Form("#pi (from D^{0}) impact parameter%s, %s, %.1f < p_{T} < %.1f GeV/c; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data(), ptLow, ptHigh), 1000, -1000., 1000.);
434  hImpParD0Pi->Sumw2();
435  hImpParD0Pi->SetLineColor(4);
436  hImpParD0Pi->SetMarkerColor(4);
437  hImpParD0Pi->SetMarkerStyle(20);
438  hImpParD0Pi->SetMarkerSize(0.6);
439  list->Add(hImpParD0Pi);
440 
441  //ImpParD0Pi2
442  TH1D* hImpParD0Pi2 = new TH1D(Form("ImpParD0Pi2%s_%d", regions[r].Data(), i), Form("#pi (from D^{0}) impact parameter (no recalc prim vtx)%s, %s, %.1f < p_{T} < %.1f GeV/c; Impact parameter [#mum]; Entries", regionTitle.Data(), listName.Data(), ptLow, ptHigh), 1000, -1000., 1000.);
443  hImpParD0Pi2->Sumw2();
444  hImpParD0Pi2->SetLineColor(4);
445  hImpParD0Pi2->SetMarkerColor(4);
446  hImpParD0Pi2->SetMarkerStyle(20);
447  hImpParD0Pi2->SetMarkerSize(0.6);
448  list->Add(hImpParD0Pi2);
449  }
450  }
451 
452  if (list == fListSignalPrompt || list == fListSignalFromB || list == fListBackground) {
453  for (Int_t i=0;i<fNPtBins;i++) {
454  TH1D* hProdd0d0All = new TH1D(Form("Prodd0d0All_%d", i), Form("d0d0, no cut, %s; d0d0 [cm^{2}]; Entries", listName.Data()), 2000, -0.0002, 0.0002);
455  hProdd0d0All->Sumw2();
456  hProdd0d0All->SetLineColor(4);
457  hProdd0d0All->SetMarkerColor(4);
458  hProdd0d0All->SetMarkerStyle(20);
459  hProdd0d0All->SetMarkerSize(0.6);
460  list->Add(hProdd0d0All);
461 
462  TH1D* hCosPntAngleAll = new TH1D(Form("CosPntAngleAll_%d", i), Form("cos #theta_{point}, no cut, %s; cos #theta_{point}; Entries", listName.Data()), 2000, -1, 1);
463  hCosPntAngleAll->Sumw2();
464  hCosPntAngleAll->SetLineColor(4);
465  hCosPntAngleAll->SetMarkerColor(4);
466  hCosPntAngleAll->SetMarkerStyle(20);
467  hCosPntAngleAll->SetMarkerSize(0.6);
468  list->Add(hCosPntAngleAll);
469 
470  TH1D* hNormDecLenAll = new TH1D(Form("NormDecLenAll_%d", i), Form("Normalized decay length, no cut, %s; Normalized decay length; Entries", listName.Data()), 2000, 0, 20);
471  hNormDecLenAll->Sumw2();
472  hNormDecLenAll->SetLineColor(4);
473  hNormDecLenAll->SetMarkerColor(4);
474  hNormDecLenAll->SetMarkerStyle(20);
475  hNormDecLenAll->SetMarkerSize(0.6);
476  list->Add(hNormDecLenAll);
477  }
478  }
479 
480  // Mother and grandmother particle IDs for D* daughters, background, peak region
481  if (list == fListBackground) {
482  Int_t sparseBins[3] = {1000, 1000, 1000}; Double_t sparseLow[3] = {-0.5, -0.5, -0.5}; Double_t sparseHigh[3] = {999.5, 999.5, 999.5};
483  THnSparseF* hPDG = new THnSparseF("PDGMotherPeak", "PDG mother, peak region, background", 3, sparseBins, sparseLow, sparseHigh);
484  list->Add(hPDG);
485  }
486 
487  // Store true impact parameter distribution
488  if (list == fListSignalFromB) {
489  for (Int_t i=0;i<fNPtBins;i++) {
490  Float_t ptLow = ptLimits[i];
491  Float_t ptHigh = ptLimits[(i+1)];
492 
493  TH1D* hImpParTrue = new TH1D(Form("ImpParTrue_%d", i), Form("True D* impact parameter, %s, %.1f < p_{T} < %.1f GeV/c; Impact parameter [#mum]; Entries", listName.Data(), ptLow, ptHigh), 1000, -1000., 1000.);
494  hImpParTrue->Sumw2();
495  hImpParTrue->SetLineColor(4);
496  hImpParTrue->SetMarkerColor(4);
497  hImpParTrue->SetMarkerStyle(20);
498  hImpParTrue->SetMarkerSize(0.6);
499  list->Add(hImpParTrue);
500  }
501  }
502 
503  char mergeMode(102); // f
504 
505  TParameter<bool> *pSingleSideband = new TParameter<bool>("SingleSideband", fSingleSideband, mergeMode);
506  list->Add(pSingleSideband);
507 
508  for (Int_t i=0;i<fNPtBins;i++) {
509  TParameter<double> *pPeakCut = new TParameter<double>(Form("PeakCut_%d", i), fPeakCut[i], mergeMode);
510  list->Add(pPeakCut);
511 
512  TParameter<double> *pSidebandCut = new TParameter<double>(Form("SidebandCut_%d", i), fSidebandCut[i], mergeMode);
513  list->Add(pSidebandCut);
514 
515  TParameter<double> *pSidebandWindow = new TParameter<double>(Form("SidebandWindow_%d", i), fSidebandWindow[i], mergeMode);
516  list->Add(pSidebandWindow);
517  }
518 }
519 
521 { // Execute analysis for current event
522  if (!fInputEvent) {
523  Error("UserExec", "NO EVENT FOUND!");
524  return;
525  }
526 
527  AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
528  TClonesArray *arrayCand = 0;
529 
530  fNEvents->Fill(0);
531 
532  if (!aodEvent && AODEvent() && IsStandardAOD()) {
533  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
534  AliAODHandler* aodHandler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
535  if(aodHandler->GetExtensions()) {
536  AliAODExtension *ext = (AliAODExtension*) aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
537  AliAODEvent *aodFromExt = ext->GetAOD();
538  arrayCand = (TClonesArray*) aodFromExt->GetList()->FindObject("Dstar");
539  }
540  }
541  else {
542  arrayCand = (TClonesArray*) aodEvent->GetList()->FindObject("Dstar");
543  }
544 
545  Int_t pdgDgDStartoD0pi[2] = {421, 211};
546  Int_t pdgDgD0toKpi[2] = {321, 211};
547 
548  TClonesArray *mcArray = 0;
549  AliAODMCHeader *mcHeader = 0;
550  if (fReadMC) {
551  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
552  if (!mcArray) {
553  AliError("Could not find Monte-Carlo in AOD");
554  return;
555  }
556  // Load MC header
557  mcHeader = (AliAODMCHeader*) aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
558  if (!mcHeader) {
559  printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
560  return;
561  }
562  }
563 
564  if (!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField()) < 0.001) {
565  return;
566  }
567  fNEvents->Fill(1);
568 
569  fCounter->StoreEvent(aodEvent, fCuts, fReadMC);
570 
571  fMagneticField = aodEvent->GetMagneticField();
572 
573  if (!fCuts->IsEventSelected(aodEvent)) {
574  return;
575  }
576  fNEvents->Fill(2);
577 
578  AliAODVertex *vtx1 = (AliAODVertex*) aodEvent->GetPrimaryVertex();
579  if (!vtx1 || vtx1->GetNContributors() < 1) {
580  return;
581  }
582  fNEvents->Fill(3);
583 
584  if (!arrayCand) {
585  AliInfo("Could not find array of HF vertices, skipping the event");
586  return;
587  }
588  else {
589  AliDebug(2, Form("Found %d vertices", arrayCand->GetEntriesFast()));
590  }
591 
592  Int_t ptBin;
593  Double_t pt;
594 
595  Int_t nSelectedProd = 0;
596  Int_t nSelectedAna = 0;
597 
598  Int_t nCand = arrayCand->GetEntriesFast();
599  for (Int_t iCand = 0; iCand<nCand; iCand++) {
600  // AliAODVertex *origVtx;
601 
602  fIsSignal = kFALSE;
603  fIsSignalPrompt = kFALSE;
604  fIsSignalFromB = kFALSE;
605  fIsBackground = kFALSE;
606 
607  AliAODRecoCascadeHF* cand = (AliAODRecoCascadeHF*) arrayCand->At(iCand);
608  AliAODRecoDecayHF2Prong *candD0 = cand->Get2Prong();
609  AliAODTrack *candSoftPi = cand->GetBachelor();
610 
611  ptBin = fCuts->PtBin(cand->Pt());
612 
613  fNEvents->Fill(4);
614 
615  if (ptBin < 0 || ptBin >= fNPtBins) {
616  continue;
617  }
618 
619  // Save some calculation effort by putting D* invariant mass cut
620  Double_t deltaInvMass = cand->DeltaInvMass();
621  if (TMath::Abs(deltaInvMass-fPDGMDStarD0) > 0.03) {
622  continue;
623  }
624 
625  Int_t isTkSelected = fCuts->IsSelected(cand, AliRDHFCuts::kTracks);
626  if (!isTkSelected) {
627  continue;
628  }
629  fNEvents->Fill(5);
630 
631  if (!fCuts->IsInFiducialAcceptance(cand->Pt(), cand->YDstar())) {
632  continue;
633  }
634  fNEvents->Fill(6);
635 
636  nSelectedProd++;
637  nSelectedAna++;
638 
639  // Match candidate to MC
640  AliAODMCParticle *partDSt = 0;
641  if (fReadMC) {
642  Int_t mcLabel = cand->MatchToMC(413, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, mcArray);
643  if (mcLabel >= 0) {
644  fIsSignal = kTRUE;
645 
646  fNEvents->Fill(10);
647 
648  partDSt = (AliAODMCParticle*) mcArray->At(mcLabel);
649 
650  if (fSkipHijing && IsFromHijing(mcArray, partDSt)) {
651  fNEvents->Fill(11);
652  continue;
653  }
654 
655  if (IsFromB(mcArray, partDSt)) {
656  fIsSignalFromB = kTRUE;
657  }
658  else {
659  fIsSignalPrompt = kTRUE;
660  }
661  }
662  else {
663  fIsBackground = kTRUE;
664  }
665  }
666 
667  // Fill histograms for topological variables before applying cuts
668  if (fReadMC) {
669  Double_t prodd0d0 = candD0->Prodd0d0();
670  Double_t cosPntAngle = candD0->CosPointingAngle();
671  Double_t normDecLen = candD0->NormalizedDecayLengthXY()*(candD0->P()/candD0->Pt());
672 
673  if (fIsSignalPrompt) {
674  ((TH1D*) fListSignalPrompt->FindObject(Form("Prodd0d0All_%d", ptBin)))->Fill(prodd0d0);
675  ((TH1D*) fListSignalPrompt->FindObject(Form("CosPntAngleAll_%d", ptBin)))->Fill(cosPntAngle);
676  ((TH1D*) fListSignalPrompt->FindObject(Form("NormDecLenAll_%d", ptBin)))->Fill(normDecLen);
677  }
678 
679  if (fIsSignalFromB) {
680  ((TH1D*) fListSignalFromB->FindObject(Form("Prodd0d0All_%d", ptBin)))->Fill(prodd0d0);
681  ((TH1D*) fListSignalFromB->FindObject(Form("CosPntAngleAll_%d", ptBin)))->Fill(cosPntAngle);
682  ((TH1D*) fListSignalFromB->FindObject(Form("NormDecLenAll_%d", ptBin)))->Fill(normDecLen);
683  }
684 
685  if (fIsBackground) {
686  ((TH1D*) fListBackground->FindObject(Form("Prodd0d0All_%d", ptBin)))->Fill(prodd0d0);
687  ((TH1D*) fListBackground->FindObject(Form("CosPntAngleAll_%d", ptBin)))->Fill(cosPntAngle);
688  ((TH1D*) fListBackground->FindObject(Form("NormDecLenAll_%d", ptBin)))->Fill(normDecLen);
689  }
690  }
691 
692  Int_t isSelected = fCuts->IsSelected(cand, AliRDHFCuts::kCandidate, aodEvent);
693  if (!isSelected) {
694  continue;
695  }
696  fNEvents->Fill(7);
697 
698  if(!aodEvent->GetPrimaryVertex()){
699  delete fNewPrimVtx; fNewPrimVtx = 0;
700  continue;
701  }
702  fNewPrimVtx = new AliAODVertex(*aodEvent->GetPrimaryVertex());
703 
704  //D* vertexing
706  if (!fDStarVtx) {
707  delete fNewPrimVtx; fNewPrimVtx = 0;
708  continue;
709  }
710 
711  fNEvents->Fill(8);
712 
713  Double_t d0, d0Err;
714  if (CalculateImpactParameter(candSoftPi, d0, d0Err)) {
715  d0 = TMath::Abs(d0*1.e4);
716  if (d0 < fImpParCut) {
717  delete fNewPrimVtx; fNewPrimVtx = 0;
718  delete fDStarVtx; fDStarVtx = 0;
719  continue;
720  }
721  }
722  if (CalculateImpactParameter((AliAODTrack*) candD0->GetDaughter(0), d0, d0Err)) {
723  d0 = TMath::Abs(d0*1.e4);
724  if (d0 < fImpParCut) {
725  delete fNewPrimVtx; fNewPrimVtx = 0;
726  delete fDStarVtx; fDStarVtx = 0;
727  continue;
728  }
729  }
730  if (CalculateImpactParameter((AliAODTrack*) candD0->GetDaughter(1), d0, d0Err)) {
731  d0 = TMath::Abs(d0*1.e4);
732  if (d0 < fImpParCut) {
733  delete fNewPrimVtx; fNewPrimVtx = 0;
734  delete fDStarVtx; fDStarVtx = 0;
735  continue;
736  }
737  }
738 
739  fNEvents->Fill(9);
740 
741  fTrueImpParForTree = fReadMC && fIsSignal ? CalculateTrueImpactParameterDStar(mcHeader, mcArray, cand) : -9998.;
742 
743  CheckInvMassDStar(cand);
744  FillHistograms(cand);
745 
746  if (fReadMC && fIsBackground && fIsPeak) {
747  AliAODTrack *daughters[3]; //0=soft pion, 1=D0 pion, 2=D0 kaon
748  daughters[0] = candSoftPi;
749 
750  if (((AliAODTrack*) candD0->GetDaughter(0))->Charge()*cand->Charge() > 0) {
751  // D* charge and first daughter charge are equal in sign <=> first daughter = pion
752  daughters[1] = (AliAODTrack*) candD0->GetDaughter(0);
753  daughters[2] = (AliAODTrack*) candD0->GetDaughter(1);
754  }
755  else {
756  daughters[1] = (AliAODTrack*) candD0->GetDaughter(1);
757  daughters[2] = (AliAODTrack*) candD0->GetDaughter(0);
758  }
759 
760  Double_t PDGCodes[3] = {0., 0., 0.};
761 
762  for (Int_t iD=0; iD<3; iD++) {
763  AliAODTrack *daughter = daughters[iD];
764 
765  Int_t daughterLabel = daughter->GetLabel();
766  if (daughterLabel >= 0) {
767  AliAODMCParticle *daughterMC = (AliAODMCParticle*) mcArray->At(daughterLabel);
768  Int_t motherLabel = daughterMC->GetMother();
769  if (motherLabel >= 0) {
770  AliAODMCParticle *motherMC = (AliAODMCParticle*) mcArray->At(motherLabel);
771  Int_t motherPDG = TMath::Abs(motherMC->GetPdgCode());
772 
773  if (iD == 0) {
774  //Soft pion: store ID of mother
775  PDGCodes[iD] = motherPDG;
776  }
777  else {
778  //D0 daughters: look at grandmother
779 
780  Int_t grandmotherLabel = motherMC->GetMother();
781  if (grandmotherLabel >= 0) {
782  AliAODMCParticle *grandmotherMC = (AliAODMCParticle*) mcArray->At(grandmotherLabel);
783  Int_t grandmotherPDG = TMath::Abs(grandmotherMC->GetPdgCode());
784  PDGCodes[iD] = grandmotherPDG;
785  }
786  }
787  }
788  }
789  }
790 
791  ((THnSparseF*) fListBackground->FindObject("PDGMotherPeak"))->Fill(PDGCodes);
792  }
793 
794  if (fReadMC && fIsSignalFromB) {
796  }
797 
798 
799  /*cand->UnsetOwnPrimaryVtx();
800  if (origVtx) {
801  cand->SetOwnPrimaryVtx(origVtx);
802  }*/
803 
804  delete fNewPrimVtx; fNewPrimVtx = 0;
805  delete fDStarVtx; fDStarVtx = 0;
806  //delete origVtx; origVtx = 0;
807  }
808 
809  fCounter->StoreCandidates(aodEvent, nSelectedProd, kTRUE);
810  fCounter->StoreCandidates(aodEvent, nSelectedAna, kFALSE);
811 
812  PostData(1, fNEvents);
813  PostData(2, fListCandidate);
814  PostData(3, fListSignal);
815  PostData(4, fListSignalPrompt);
816  PostData(5, fListSignalFromB);
817  PostData(6, fListBackground);
818  PostData(8, fCounter);
819  PostData(9, fTreeCandidate);
820 }
821 
823 { // Determine primary vertex without D* daughters
824  AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
825  if (!vtxAOD) {
826  return 0;
827  }
828 
829  TString title = vtxAOD->GetTitle();
830  if (!title.Contains("VertexerTracks")) {
831  return 0;
832  }
833 
834  AliVertexerTracks *vertexer = new AliVertexerTracks(aod->GetMagneticField());
835 
836  vertexer->SetITSMode();
837  vertexer->SetMinClusters(3);
838  vertexer->SetConstraintOff();
839 
840  if (title.Contains("WithConstraint")) {
841  Float_t diamondcovxy[3];
842  aod->GetDiamondCovXY(diamondcovxy);
843  Double_t pos[3] = {aod->GetDiamondX(), aod->GetDiamondY(), 0.};
844  Double_t cov[6] = {diamondcovxy[0], diamondcovxy[1], diamondcovxy[2], 0., 0., 10.*10.};
845  AliESDVertex *diamond = new AliESDVertex(pos, cov, 1., 1);
846  vertexer->SetVtxStart(diamond);
847  delete diamond; diamond=NULL;
848  }
849 
850  Int_t skipped[10]; for(Int_t i=0; i<10; i++) { skipped[i] = -1; }
851  Int_t nTrksToSkip = 0, id;
852  AliAODTrack *t = 0;
853 
854  for (Int_t i=0; i<3; i++) {
855  if (i == 0) {
856  t = cand->GetBachelor();
857  }
858  else {
859  AliAODRecoDecayHF2Prong *twoProng = cand->Get2Prong();
860  t = (AliAODTrack*) twoProng->GetDaughter(i-1);
861  }
862 
863  id = (Int_t) t->GetID();
864  if (id < 0) {
865  continue;
866  }
867  skipped[nTrksToSkip++] = id;
868  }
869 
870  vertexer->SetSkipTracks(nTrksToSkip, skipped);
871  AliESDVertex *vtxESDNew = vertexer->FindPrimaryVertex(aod);
872 
873  delete vertexer; vertexer = NULL;
874 
875  if (!vtxESDNew) {
876  return 0;
877  }
878 
879  if (vtxESDNew->GetNContributors() <= 0) {
880  delete vtxESDNew; vtxESDNew = NULL;
881  return 0;
882  }
883 
884  // convert to AliAODVertex
885  Double_t pos[3], cov[6], chi2perNDF;
886  vtxESDNew->GetXYZ(pos); // position
887  vtxESDNew->GetCovMatrix(cov); //covariance matrix
888  chi2perNDF = vtxESDNew->GetChi2toNDF();
889  delete vtxESDNew; vtxESDNew = NULL;
890 
891  AliAODVertex *vtxAODNew = new AliAODVertex(pos, cov, chi2perNDF);
892 
893  return vtxAODNew;
894 }
895 
897 { // Determine the D* vertex
898 
899  // Double_t dcadummy, xdummy, ydummy;
900  AliAODRecoDecayHF2Prong *candD0 = cand->Get2Prong();
901  AliAODTrack *candSoftPi = cand->GetBachelor();
902  AliNeutralTrackParam *ntpD0 = new AliNeutralTrackParam(candD0);
903  // ntpD0->PropagateToDCA(fNewPrimVtx, fMagneticField, kVeryBig);
904  //AliExternalTrackParam etpSoftPi; etpSoftPi.CopyFromVTrack(candSoftPi);
905  AliESDtrack *candSoftPiESD = new AliESDtrack(candSoftPi);
906  // dcadummy = candSoftPiESD->GetDCA(ntpD0, fMagneticField, xdummy, ydummy);
907 
908  TObjArray *tracks = new TObjArray(2);
909  tracks->AddAt(candSoftPiESD, 0);
910  tracks->AddAt(ntpD0, 1);
911  // tracks->Add(&etpSoftPi);
912 
913  Double_t pos[3], cov[6];
914  fNewPrimVtx->GetXYZ(pos);
915  fNewPrimVtx->GetCovarianceMatrix(cov);
916  AliESDVertex *newPrimVtxESD = new AliESDVertex(pos, cov, 100., 100, fNewPrimVtx->GetName());
917 
918  AliVertexerTracks *vertexerTracks = new AliVertexerTracks(fMagneticField);
919  vertexerTracks->SetVtxStart(newPrimVtxESD);
920  AliESDVertex *vertexESD = (AliESDVertex*) vertexerTracks->VertexForSelectedESDTracks(tracks);
921  if (!vertexESD) {
922  return 0;
923  }
924 
925  if (vertexESD->GetNContributors() != 2) {
926  delete vertexESD; vertexESD = 0;
927  return 0;
928  }
929 
930  Double_t vertRadius2 = vertexESD->GetX()*vertexESD->GetX() + vertexESD->GetY()*vertexESD->GetY();
931  if (vertRadius2 > 8.) {
932  // vertex outside beam pipe, reject candidate to avoid propagation through material
933  delete vertexESD; vertexESD = 0;
934  return 0;
935  }
936 
937  Double_t pos2[3], cov2[6], chi2perNDF;
938  vertexESD->GetXYZ(pos2); // position
939  vertexESD->GetCovMatrix(cov2); //covariance matrix
940  chi2perNDF = vertexESD->GetChi2toNDF();
941 
942  delete vertexESD; vertexESD = 0;
943 
944  AliAODVertex *vertexAOD = new AliAODVertex(pos2, cov2, chi2perNDF, 0x0, -1, AliAODVertex::kUndef, 2);
945 
946  return vertexAOD;
947 }
948 
950 { // Check if D* candidate falls within peak or sideband region
951  Double_t invMassDiff = cand->DeltaInvMass()-fPDGMDStarD0;
952  Double_t invMassDiffAbs = TMath::Abs(invMassDiff);
953 
954  fIsPeak = kFALSE;
955  fIsSideband = kFALSE;
956 
957  Int_t ptBin = fCuts->PtBin(cand->Pt());
958  if (ptBin < 0) {
959  return;
960  }
961 
962  Double_t peakCut = fPeakCut[ptBin];
963  Double_t sidebandCut = fSidebandCut[ptBin];
964  Double_t sidebandWindow = fSidebandWindow[ptBin];
965 
966  if (invMassDiffAbs < peakCut) {
967  fIsPeak = kTRUE;
968  return;
969  }
970 
971  if (!fSingleSideband) {
972  invMassDiff = invMassDiffAbs;
973  }
974 
975  if (invMassDiff > sidebandCut && invMassDiff <= sidebandCut+sidebandWindow) {
976  fIsSideband = kTRUE;
977  }
978 }
979 
980 Bool_t AliAnalysisTaskSEDStarCharmFraction::IsFromB(TClonesArray *arrayMC, const AliAODMCParticle *mcPartCandidate)
981 { // Check if the MC particle comes from a B
982  Int_t mother = mcPartCandidate->GetMother();
983  while (mother > 0) {
984  AliAODMCParticle *mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
985  if (mcGranma) {
986  Int_t pdgGranma = mcGranma->GetPdgCode();
987  Int_t abspdgGranma = TMath::Abs(pdgGranma);
988  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
989  return kTRUE;
990  }
991  mother = mcGranma->GetMother();
992  }
993  else{
994  AliError("Failed casting the mother particle!");
995  break;
996  }
997  }
998 
999  return kFALSE;
1000 }
1001 
1002 Bool_t AliAnalysisTaskSEDStarCharmFraction::IsFromHijing(TClonesArray *arrayMC, const AliAODMCParticle *mcPartCandidate)
1003 { // Check if the MC particle is from Hijing
1004  Int_t mother = mcPartCandidate->GetMother();
1005  while (mother > 0) {
1006  AliAODMCParticle *mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1007  if (mcGranma) {
1008  Int_t pdgGranma = mcGranma->GetPdgCode();
1009  Int_t abspdgGranma = TMath::Abs(pdgGranma);
1010  if (abspdgGranma == 4 || abspdgGranma == 5) {
1011  return kFALSE;
1012  }
1013  mother = mcGranma->GetMother();
1014  }
1015  else{
1016  AliError("Failed casting the mother particle!");
1017  break;
1018  }
1019  }
1020 
1021  return kTRUE;
1022 }
1023 
1025 { // Calculate impact parameter for a track
1026  AliExternalTrackParam etp;
1027  etp.CopyFromVTrack(track);
1028  d0 = 0; d0Err = 0;
1029  Double_t dz[2], covdz[3];
1030  if (etp.PropagateToDCA(fNewPrimVtx, fMagneticField, 3., dz, covdz)) {
1031  d0 = dz[0];
1032  d0Err = TMath::Sqrt(covdz[0]);
1033  return kTRUE;
1034  }
1035 
1036  return kFALSE;
1037 }
1038 
1040 { // Calculate impact parameter of the D*
1041  // Calculation from AliAODRecoDecay::ImpParXY()
1042 
1043  Double_t p[3] = {cand->Px(), cand->Py(), cand->Pz()};
1044  Double_t ptsq = p[0]*p[0] + p[1]*p[1];
1045 
1046  Double_t prim[3] = {fNewPrimVtx->GetX(), fNewPrimVtx->GetY(), fNewPrimVtx->GetZ()};
1047  Double_t sec[3] = {fDStarVtx->GetX(), fDStarVtx->GetY(), fDStarVtx->GetZ()};
1048 
1049  Double_t k = - (sec[0]-prim[0])*p[0] - (sec[1]-prim[1])*p[1];
1050  k /= ptsq;
1051  Double_t dx = sec[0] - prim[0] + k*p[0];
1052  Double_t dy = sec[1] - prim[1] + k*p[1];
1053  Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
1054 
1055  TVector3 mom(p[0], p[1], p[2]);
1056  TVector3 fline(sec[0]-prim[0], sec[1]-prim[1], sec[2]-prim[2]);
1057  TVector3 cross = mom.Cross(fline);
1058 
1059  return (cross.Z() > 0. ? absImpPar : -absImpPar);
1060 }
1061 
1063 { // Calculate true impact parameter of the D*
1064  AliAODRecoDecayHF2Prong *candD0 = cand->Get2Prong();
1065  AliAODTrack *tr1 = (AliAODTrack*) candD0->GetDaughter(0);
1066  AliAODTrack *tr2 = (AliAODTrack*) candD0->GetDaughter(1);
1067  AliAODTrack *tr3 = (AliAODTrack*) cand->GetBachelor();
1068 
1069  Int_t label1 = tr1->GetLabel();
1070  Int_t label2 = tr2->GetLabel();
1071  Int_t label3 = tr3->GetLabel();
1072 
1073  if (label1 < 0 || label2 < 0 || label3 < 0) {
1074  return -9999.;
1075  }
1076 
1077  AliAODMCParticle *part1 = (AliAODMCParticle*) arrayMC->At(label1);
1078  AliAODMCParticle *part2 = (AliAODMCParticle*) arrayMC->At(label2);
1079  AliAODMCParticle *part3 = (AliAODMCParticle*) arrayMC->At(label3);
1080 
1081  if (part1 == 0 || part2 == 0 || part3 == 0) {
1082  return -9999.;
1083  }
1084 
1085  Double_t PxTotal = part1->Px()+part2->Px()+part3->Px();
1086  Double_t PyTotal = part1->Py()+part2->Py()+part3->Py();
1087  Double_t PzTotal = part1->Pz()+part2->Pz()+part3->Pz();
1088 
1089 
1090  // True primary+secondary vertices and momentum
1091  Double_t prim[3] = {headerMC->GetVtxX(), headerMC->GetVtxY(), headerMC->GetVtxZ()};
1092  Double_t sec[3] = {part3->Xv(), part3->Yv(), part3->Zv()}; // Put soft pion production vertex as secondary vertex
1093  Double_t p[3] = {PxTotal, PyTotal, PzTotal};
1094 
1095 
1096  Double_t ptsq = p[0]*p[0] + p[1]*p[1];
1097  Double_t k = - (sec[0]-prim[0])*p[0] - (sec[1]-prim[1])*p[1];
1098  k /= ptsq;
1099  Double_t dx = sec[0] - prim[0] + k*p[0];
1100  Double_t dy = sec[1] - prim[1] + k*p[1];
1101  Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy)*1.e4;
1102 
1103  TVector3 mom(p[0], p[1], p[2]);
1104  TVector3 fline(sec[0]-prim[0], sec[1]-prim[1], sec[2]-prim[2]);
1105  TVector3 cross = mom.Cross(fline);
1106 
1107  return cross.Z() > 0. ? absImpPar : -absImpPar;
1108 }
1109 
1111 { // Fill histograms for a D* candidate
1112  Double_t Pt = cand->Pt();
1113  Int_t ptBin = fCuts->PtBin(Pt);
1114  Double_t deltaInvMass = cand->DeltaInvMass();
1115  Double_t impPar = CalculateImpactParameterDStar(cand)*1.e4;
1116  Double_t invMassD0 = cand->InvMassD0();
1117 
1118  FillHistogram("DeltaInvMass", deltaInvMass);
1119  if (ptBin >= 0) {
1120  FillHistogram(Form("DeltaInvMass_%d", ptBin), deltaInvMass);
1121  }
1122 
1123  FillHistogram("InvMassD0", invMassD0);
1124  if (ptBin >= 0) {
1125  FillHistogram(Form("InvMassD0_%d", ptBin), invMassD0);
1126  }
1127 
1128  FillRegionHistogram("ImpPar$r", impPar);
1129  if (ptBin >= 0) {
1130  FillRegionHistogram(Form("ImpPar$r_%d", ptBin), impPar);
1131  }
1132 
1133  // Store candidate in tree
1134  fPtForTree = Pt;
1135  fInvMassForTree = deltaInvMass;
1136  fImpParForTree = impPar;
1137  // fTrueImpParForTree was already updated before this
1138  fTypeForTree = fIsPeak ? 1 : (fIsSideband ? (deltaInvMass < fPDGMDStarD0 ? 2 : 3) : 0);
1139  fSignalTypeForTree = fReadMC ? (fIsSignalPrompt ? 1 : (fIsSignalFromB ? 2 : 0)) : -1;
1140  fTreeCandidate->Fill();
1141 
1142  AliAODTrack *softPi = cand->GetBachelor();
1143  Double_t d0, d0Err;
1144  if (CalculateImpactParameter(softPi, d0, d0Err)) {
1145  FillRegionHistogram("ImpParSoftPi$r", d0*1.e4);
1146  if (ptBin >= 0) {
1147  FillRegionHistogram(Form("ImpParSoftPi$r_%d", ptBin), d0*1.e4);
1148  }
1149  }
1150 
1151  AliAODRecoDecayHF2Prong *candD0 = cand->Get2Prong();
1152  AliAODTrack *pi; AliAODTrack *K; Int_t piId, KId;
1153  if (((AliAODTrack*) candD0->GetDaughter(0))->Charge()*cand->Charge() > 0) {
1154  // D* and first daughter charge are equal in sign <=> first daughter = pion
1155  piId = 0; KId = 1;
1156  }
1157  else {
1158  piId = 1; KId = 0;
1159  }
1160 
1161  pi = (AliAODTrack*) candD0->GetDaughter(piId);
1162  K = (AliAODTrack*) candD0->GetDaughter(KId);
1163 
1164  //pi (from D0) impact parameter
1165  if (CalculateImpactParameter(pi, d0, d0Err)) {
1166  FillRegionHistogram("ImpParD0Pi$r", d0*1.e4);
1167  if (ptBin >= 0) {
1168  FillRegionHistogram(Form("ImpParD0Pi$r_%d", ptBin), d0*1.e4);
1169  }
1170  }
1171 
1172  //pi (from D0) impact parameter without recalculation of primary vertex
1173  d0 = candD0->Getd0Prong(piId);
1174  FillRegionHistogram("ImpParD0Pi2$r", d0*1.e4);
1175  if (ptBin >= 0) {
1176  FillRegionHistogram(Form("ImpParD0Pi2$r_%d", ptBin), d0*1.e4);
1177  }
1178 
1179  //K (from D0) impact parameter
1180  if (CalculateImpactParameter(K, d0, d0Err)) {
1181  FillRegionHistogram("ImpParD0K$r", d0*1.e4);
1182  if (ptBin >= 0) {
1183  FillRegionHistogram(Form("ImpParD0K$r_%d", ptBin), d0*1.e4);
1184  }
1185  }
1186 
1187  //K (from D0) impact parameter without recalculation of primary vertex
1188  d0 = candD0->Getd0Prong(KId);
1189  FillRegionHistogram("ImpParD0K2$r", d0*1.e4);
1190  if (ptBin >= 0) {
1191  FillRegionHistogram(Form("ImpParD0K2$r_%d", ptBin), d0*1.e4);
1192  }
1193 
1194  AliAODVertex *origVtx = 0;
1195  if (candD0->GetOwnPrimaryVtx()) {
1196  origVtx = new AliAODVertex(*candD0->GetOwnPrimaryVtx());
1197  }
1198  candD0->SetOwnPrimaryVtx(fNewPrimVtx);
1199 
1200  Double_t impParD0 = candD0->ImpParXY()*1.e4;
1201  FillRegionHistogram("ImpParD0$r", impParD0);
1202  if (ptBin >= 0) {
1203  FillRegionHistogram(Form("ImpParD0$r_%d", ptBin), impParD0);
1204  }
1205 
1206  candD0->UnsetOwnPrimaryVtx();
1207  if (origVtx) {
1208  candD0->SetOwnPrimaryVtx(origVtx);
1209  }
1210  delete origVtx;
1211 }
1212 
1214 { // Fill a specific histogram in multiple lists
1215  ((TH1D*) fListCandidate->FindObject(name))->Fill(value);
1216 
1217  if (fIsSignal) {
1218  ((TH1D*) fListSignal->FindObject(name))->Fill(value);
1219  }
1220 
1221  if (fIsSignalPrompt) {
1222  ((TH1D*) fListSignalPrompt->FindObject(name))->Fill(value);
1223  }
1224 
1225  if (fIsSignalFromB) {
1226  ((TH1D*) fListSignalFromB->FindObject(name))->Fill(value);
1227  }
1228 
1229  if (fIsBackground) {
1230  ((TH1D*) fListBackground->FindObject(name))->Fill(value);
1231  }
1232 }
1233 
1235 { // Fill a specific histogram in multiple lists, in the approprate regions (all, peak region, sideband region)
1236  FillHistogram(TString(name).ReplaceAll("$r", "").Data(), value);
1237 
1238  if (fIsPeak) {
1239  FillHistogram(TString(name).ReplaceAll("$r", "Peak").Data(), value);
1240  }
1241 
1242  if (fIsSideband) {
1243  FillHistogram(TString(name).ReplaceAll("$r", "Sideband").Data(), value);
1244  }
1245 }
1246 
1248 { // Fill histogram with true impact parameter distribution for D from B
1249  Double_t Pt = cand->Pt();
1250  Int_t ptBin = fCuts->PtBin(Pt);
1251 
1252  if (ptBin < 0) {
1253  return;
1254  }
1255 
1256  if (fTrueImpParForTree == -9999.) {
1257  return;
1258  }
1259 
1260  TH1D* hImpParTrue = (TH1D*) fListSignalFromB->FindObject(Form("ImpParTrue_%d", ptBin));
1261  hImpParTrue->Fill(fTrueImpParForTree);
1262 }
1263 
1265 { // Terminate
1266 }
Double_t NormalizedDecayLengthXY() const
Bool_t IsFromHijing(TClonesArray *arrayMC, const AliAODMCParticle *mcPartCandidate)
double Double_t
Definition: External.C:58
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
const char * title
Definition: MakeQAPdf.C:27
Double_t DeltaInvMass() const
void FillRegionHistogram(const char *name, Double_t value)
Bool_t CalculateImpactParameter(AliAODTrack *track, Double_t &d0, Double_t &d0Err)
Double_t YDstar() const
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
Double_t ImpParXY() const
void FillHistogram(const char *name, Double_t value)
Double_t CalculateImpactParameterDStar(AliAODRecoCascadeHF *cand)
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Definition: External.C:212
AliAODTrack * GetBachelor() const
AliAODVertex * GetOwnPrimaryVtx() const
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod)
AliAODVertex * RemoveDaughtersFromPrimaryVtx(AliAODEvent *aod, AliAODRecoCascadeHF *cand)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Bool_t IsEventSelected(AliVEvent *event)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
AliAODVertex * ReconstructDStarVtx(AliAODRecoCascadeHF *cand)
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
Float_t * GetPtBinLimits() const
Definition: AliRDHFCuts.h:248
Double_t CalculateTrueImpactParameterDStar(AliAODMCHeader *headerMC, TClonesArray *arrayMC, AliAODRecoCascadeHF *cand)
Double_t InvMassD0() const
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:249
Bool_t IsFromB(TClonesArray *arrayMC, const AliAODMCParticle *mcPartCandidate)
const Double_t pi
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
AliAODRecoDecayHF2Prong * Get2Prong() const
Int_t PtBin(Double_t pt) const