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