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