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