AliPhysics  96f6795 (96f6795)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskSEDstoK0sK.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2010, 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 /* $Id$ */
17 
31 
32 
33 #include <TClonesArray.h>
34 #include <TDatabasePDG.h>
35 #include <TLorentzVector.h>
36 #include <TVector3.h>
37 
38 #include "AliAnalysisManager.h"
39 #include "AliAnalysisTaskSE.h"
41 #include "AliAnalysisVertexingHF.h"
42 #include "AliAODEvent.h"
43 #include "AliAODHandler.h"
44 #include "AliAODMCHeader.h"
45 #include "AliAODMCParticle.h"
46 #include "AliAODPidHF.h"
47 #include "AliAODRecoCascadeHF.h"
48 #include "AliESDtrack.h"
49 #include "AliESDVertex.h"
51 #include "AliNeutralTrackParam.h"
52 #include "AliRDHFCuts.h"
53 #include "AliRDHFCutsDstoK0sK.h"
54 #include "AliVertexerTracks.h"
55 #include "AliVTrack.h"
56 
57 
61 
62 
63 
64 
65 
66 //__________________________________________________________________________
69  , fAnalysisCuts(0)
70  , fCounter(0)
71  , fOutputSele(0)
72  , fOutputCand(0)
73  , fOutputPID(0)
74  , fOutputMC(0)
75  , fOutputNtuple(0)
76  , fHisNEvents(0)
77  , fHisRapidity(0)
78  , fHisRapiditySel(0)
79  , fHisPidTPCKaonVsPt(0)
80  , fHisPidTOFKaonVsPt(0)
81  , fHisPidTPCTOFKaon(0)
82  , fHisPidTPCTOFKaonSel(0)
83  , fHisPidTPCKaonVsPtSel(0)
84  , fHisPidTOFKaonVsPtSel(0)
85  , fNPtBins(0)
86  , fMassRange(0.8)
87  , fMassBinSize(0.002)
88  , fReadMC(kFALSE)
89  , fUseSelectionBit(kFALSE)
90  , fFillNtuple(kFALSE)
91  , fAODProtection(0)
92 {
96 
97  for (Int_t ihis=0; ihis<3; ihis++) {
98  fHisCentrality[ihis] = 0;
99  fHisCentralityVsMult[ihis] = 0;
100  }
101 
102  for (Int_t ihis=0; ihis<5; ihis++) {
103  fHisInvMassDs[ihis] = 0;
104  fHisInvMassDplus[ihis] = 0;
105  fHisInvMassK0s[ihis] = 0;
106  fHisPtK0s[ihis] = 0;
107  fHisPtBachelor[ihis] = 0;
108  fHisImpParK0s[ihis] = 0;
109  fHisImpParBach[ihis] = 0;
110  fHisCTauK0s[ihis] = 0;
111  fHisCosPointingDs[ihis] = 0;
112  fHisCosPointingXYDs[ihis] = 0;
113  fHisCosThetaStarK0s[ihis] = 0;
114  fHisCosThetaStarBach[ihis] = 0;
115  fHisDCAK0sBach[ihis] = 0;
116  fHisDecayLxyDs[ihis] = 0;
117  fHisNormDecayLxyDs[ihis] = 0;
118  }
119 
120  for (Int_t ihis=0; ihis<kMaxPtBins+1; ihis++) {
121  fPtLimits[ihis] = 0;
122  }
123 
124  for (Int_t ihis=0; ihis<kNTupleVars; ihis++) {
125  fCutsMinTupleVars[ihis] = 0;
126  fCutsMaxTupleVars[ihis] = 0;
127  }
128 }
129 
130 
131 
132 
133 
134 //__________________________________________________________________________
136  AliRDHFCutsDstoK0sK *cuts,
137  Bool_t readMC,
138  Bool_t fillNtuple,
139  Int_t nCutsTuple,
140  Float_t *minCutsTuple,
141  Float_t *maxCutsTuple) :
142 AliAnalysisTaskSE(name)
143  , fAnalysisCuts(cuts)
144  , fCounter(0)
145  , fOutputSele(0)
146  , fOutputCand(0)
147  , fOutputPID(0)
148  , fOutputMC(0)
149  , fOutputNtuple(0)
150  , fHisNEvents(0)
151  , fHisRapidity(0)
152  , fHisRapiditySel(0)
153  , fHisPidTPCKaonVsPt(0)
154  , fHisPidTOFKaonVsPt(0)
155  , fHisPidTPCTOFKaon(0)
156  , fHisPidTPCTOFKaonSel(0)
157  , fHisPidTPCKaonVsPtSel(0)
158  , fHisPidTOFKaonVsPtSel(0)
159  , fNPtBins(0)
160  , fMassRange(0.8)
161  , fMassBinSize(0.002)
162  , fReadMC(readMC)
163  , fUseSelectionBit(kFALSE)
164  , fFillNtuple(fillNtuple)
165  , fAODProtection(0)
166 {
170 
171  for (Int_t ihis=0; ihis<3; ihis++) {
172  fHisCentrality[ihis] = 0;
173  fHisCentralityVsMult[ihis] = 0;
174  }
175 
176  for (Int_t ihis=0; ihis<5;ihis++) {
177  fHisInvMassDs[ihis] = 0;
178  fHisInvMassDplus[ihis] = 0;
179  fHisInvMassK0s[ihis] = 0;
180  fHisPtK0s[ihis] = 0;
181  fHisPtBachelor[ihis] = 0;
182  fHisImpParK0s[ihis] = 0;
183  fHisImpParBach[ihis] = 0;
184  fHisCTauK0s[ihis] = 0;
185  fHisCosPointingDs[ihis] = 0;
186  fHisCosPointingXYDs[ihis] = 0;
187  fHisCosThetaStarK0s[ihis] = 0;
188  fHisCosThetaStarBach[ihis] = 0;
189  fHisDCAK0sBach[ihis] = 0;
190  fHisDecayLxyDs[ihis] = 0;
191  fHisNormDecayLxyDs[ihis] = 0;
192  }
193 
194  for (Int_t ihis=0; ihis<kMaxPtBins+1; ihis++) {
195  fPtLimits[ihis] = 0;
196  }
197 
198  for (Int_t ihis=0; ihis<kNTupleVars; ihis++) {
199  fCutsMinTupleVars[ihis] = 0;
200  fCutsMaxTupleVars[ihis] = 0;
201  }
202 
203 
204 
206  if (fFillNtuple)
207  SetCutsTupleVariables(nCutsTuple, minCutsTuple, maxCutsTuple);
208 
209 
210  // Output slot #1 writes into a TList container (fAnalysisCuts)
211  DefineOutput(1, TList::Class());
212 
213  // Output slot #2 writes into a AliNormalizationCounter container (fCounter)
214  DefineOutput(2, AliNormalizationCounter::Class());
215 
216  // Output slot #3 writes into a TList container (fOutputSele)
217  DefineOutput(3, TList::Class());
218 
219  if (!fFillNtuple) {
220  // Output slot #4 writes into a TList container (fOutputCand)
221  DefineOutput(4, TList::Class());
222 
223  // Output slot #5 writes into a TList container (fOutputPID)
224  DefineOutput(5, TList::Class());
225 
226  if (fReadMC) {
227  // Output slot #6 writes into a TList container (fOutputMC)
228  DefineOutput(6, TList::Class());
229  }
230  } else {
231  // Output slot #4 writes into a TNtuple container (fOutputNtuple)
232  DefineOutput(4, TNtuple::Class());
233  }
234 }
235 
236 
237 
238 
239 
240 //__________________________________________________________________________
242 {
246 
247  if (fAnalysisCuts) { delete fAnalysisCuts; fAnalysisCuts =0; }
248  if (fCounter) { delete fCounter; fCounter =0; }
249  if (fOutputSele) { delete fOutputSele; fOutputSele =0; }
250  if (fOutputCand) { delete fOutputCand; fOutputCand =0; }
251  if (fOutputPID) { delete fOutputPID; fOutputPID =0; }
252  if (fOutputMC) { delete fOutputMC; fOutputMC =0; }
253  if (fOutputNtuple) { delete fOutputNtuple; fOutputNtuple =0; }
254 }
255 
256 
257 
258 
259 
260 //__________________________________________________________________________
262 {
266  AliDebug(1, "AliAnalysisTaskSEDstoK0sK::Init()");
267 
268 
269  //---------------------------------------------------------------------------
270  // - Output slot #1: cut objects
271  //---------------------------------------------------------------------------
273  listCuts->SetName(GetOutputSlot(1)->GetContainer()->GetName());
274  PostData(1, listCuts);
275 
276  return;
277 }
278 
279 
280 
281 
282 
283 //__________________________________________________________________________
285 {
289  AliDebug(1, "AliAnalysisTaskSEDstoK0sK::UserCreateOutputObjects()");
290 
291 
292 
293  // - Define the Pt range and binning
294  Double_t ptBinsRange[fNPtBins+1];
295  for (Int_t ipt=0; ipt<fNPtBins+1; ipt++) {
296  ptBinsRange[ipt] = fPtLimits[ipt];
297  }
298 
299 
300 
301  //---------------------------------------------------------------------------
302  // - Output slot #2: counter for Normalization
303  //---------------------------------------------------------------------------
304  fCounter = new AliNormalizationCounter(GetOutputSlot(2)->GetContainer()->GetName());
305  fCounter->Init();
306 
307 
308 
309  //---------------------------------------------------------------------------
310  // - Output slot #3: various histograms of selected events
311  //---------------------------------------------------------------------------
312  fOutputSele = new TList();
313  fOutputSele->SetOwner();
314  fOutputSele->SetName(GetOutputSlot(3)->GetContainer()->GetName());
315 
316 
317  fHisNEvents = new TH1F("fHisNEvents", ";;Entries", 21, 0, 21);
318  fHisNEvents->GetXaxis()->SetBinLabel(1, "nEvents Read");
319  fHisNEvents->GetXaxis()->SetBinLabel(2, "nEvents Matched dAOD");
320  fHisNEvents->GetXaxis()->SetBinLabel(3, "nEvents Mismatched dAOD");
321  fHisNEvents->GetXaxis()->SetBinLabel(4, "Analysed events");
322  fHisNEvents->GetXaxis()->SetBinLabel(5, "Rejection: Trigger");
323  fHisNEvents->GetXaxis()->SetBinLabel(6, "Rejection: Vertex reco");
324  fHisNEvents->GetXaxis()->SetBinLabel(7, "Rejection: Pileup");
325  fHisNEvents->GetXaxis()->SetBinLabel(8, "Rejection: Centrality");
326  fHisNEvents->GetXaxis()->SetBinLabel(9, "Rejection: VtxZ fiducial acc.");
327  fHisNEvents->GetXaxis()->SetBinLabel(10, "Rejection: Physics Selection");
328  fHisNEvents->GetXaxis()->SetBinLabel(11, "Selected events");
329  fHisNEvents->GetXaxis()->SetBinLabel(12, "Read candidates");
330  fHisNEvents->GetXaxis()->SetBinLabel(13, "FillRecoCasc fails");
331  fHisNEvents->GetXaxis()->SetBinLabel(14, "Offline SV reco. fails");
332  fHisNEvents->GetXaxis()->SetBinLabel(15, "Cand with sec. vertex");
333  fHisNEvents->GetXaxis()->SetBinLabel(16, "In fiducial acceptance");
334  fHisNEvents->GetXaxis()->SetBinLabel(17, "Selected Tracks");
335  fHisNEvents->GetXaxis()->SetBinLabel(18, "Selected Candidates");
336  fHisNEvents->GetXaxis()->SetBinLabel(19, "Selected PID");
337  fHisNEvents->GetXaxis()->SetBinLabel(20, "Selected true D^{+}_{s}");
338  fHisNEvents->GetXaxis()->SetBinLabel(21, "Selected true D^{+}");
339  fOutputSele->Add(fHisNEvents);
340 
341 
342  fHisCentrality[0] = new TH1F("fHisCentrality_all", "; Centrality; Entries", 1000, 0., 100.);
343  fHisCentrality[0]->SetTitle("All read events");
344 
345  fHisCentrality[1] = (TH1F*) fHisCentrality[0]->Clone("fHisCentrality_selected");
346  fHisCentrality[1]->SetTitle("Selected events");
347 
348  fHisCentrality[2] = (TH1F*) fHisCentrality[0]->Clone("fHisCentrality_rejected");
349  fHisCentrality[2]->SetTitle("Events rejected due to centrality");
350 
351  fOutputSele->Add(fHisCentrality[0]);
352  fOutputSele->Add(fHisCentrality[1]);
353  fOutputSele->Add(fHisCentrality[2]);
354 
355 
356  fHisCentralityVsMult[0] = new TH2F("fHisCentralityVsMult_all", "; Track multiplicities; Centrality; Entries", 100, 0., 10000., 40, 0., 100.);
357  fHisCentralityVsMult[0]->SetTitle("All read events");
358 
359  fHisCentralityVsMult[1] = (TH2F*) fHisCentralityVsMult[0]->Clone("fHisCentralityVsMult_selected");
360  fHisCentralityVsMult[1]->SetTitle("Selected events");
361 
362  fHisCentralityVsMult[2] = (TH2F*) fHisCentralityVsMult[0]->Clone("fHisCentralityVsMult_rejected");
363  fHisCentralityVsMult[2]->SetTitle("Events rejected due to centrality");
364 
368 
369 
370  fHisRapidity = new TH2F("fHisRapidity_pt", "; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); #it{y}(D^{+}_{s}); Entries", 40, 0., 20., 80, -2., 2.);
371  fHisRapidity->SetTitle("Production cuts + SV reconstruction");
372 
373  fHisRapiditySel = (TH2F*) fHisRapidity->Clone("fHisRapiditySel_pt");
374  fHisRapiditySel->SetTitle("Selected candidates");
375 
378 
379 
380  fHisPidTPCKaonVsPt = new TH2F("fHisPidTPCKaonVsPt", "; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c});; Entries", fNPtBins, ptBinsRange, 200, -10., 10.);
381  fHisPidTPCKaonVsPt->GetYaxis()->SetTitle("n#it{#sigma}_{TPC}^{K}");
382  fHisPidTPCKaonVsPt->SetTitle("Selected tracks + candidates");
383 
384  fHisPidTPCKaonVsPtSel = (TH2F*) fHisPidTPCKaonVsPt->Clone("fHisPidTPCKaonVsPtSel");
385  fHisPidTPCKaonVsPtSel->SetTitle("Selected tracks + candidates + PID");
386 
389 
390 
391  fHisPidTOFKaonVsPt = (TH2F*) fHisPidTPCKaonVsPt->Clone("fHisPidTOFKaonVsPt");
392  fHisPidTOFKaonVsPt->GetYaxis()->SetTitle("n#it{#sigma}_{TOF}^{K}");
393  fHisPidTOFKaonVsPt->SetTitle("Selected tracks + candidates");
394 
395  fHisPidTOFKaonVsPtSel = (TH2F*) fHisPidTOFKaonVsPt->Clone("fHisPidTOFKaonVsPtSel");
396  fHisPidTOFKaonVsPtSel->SetTitle("Selected tracks + candidates + PID");
397 
400 
401 
402  fHisPidTPCTOFKaon = new TH2F("fHisPidTPCTOFKaon", "; n#it{#sigma}_{TPC}^{K}; n#it{#sigma}_{TOF}^{K}; Entries", 200, -10., 10., 200, -10., 10.);
403  fHisPidTPCTOFKaon->SetTitle("Selected tracks + candidates");
404 
405  fHisPidTPCTOFKaonSel = (TH2F*) fHisPidTPCTOFKaon->Clone("fHisPidTPCTOFKaonSel");
406  fHisPidTPCTOFKaonSel->SetTitle("Selected tracks + candidates + PID");
407 
410 
411 
412 
413  if (!fFillNtuple) {
414 
415  //---------------------------------------------------------------------------
416  // - Output slot #4 and #5: histograms at candidate/PID level
417  // - Output slot #6: histograms for MC informations
418  //---------------------------------------------------------------------------
419  fOutputCand = new TList();
420  fOutputCand->SetOwner();
421  fOutputCand->SetName(GetOutputSlot(4)->GetContainer()->GetName());
422 
423  fOutputPID = new TList();
424  fOutputPID->SetOwner();
425  fOutputPID->SetName(GetOutputSlot(5)->GetContainer()->GetName());
426 
427  if (fReadMC) {
428  fOutputMC = new TList();
429  fOutputMC->SetOwner();
430  fOutputMC->SetName(GetOutputSlot(6)->GetContainer()->GetName());
431  }
432 
433 
434  // - Define the invariant mass range and binning
435  Float_t massDs = TDatabasePDG::Instance()->GetParticle(431)->Mass();
436  Float_t massDplus = TDatabasePDG::Instance()->GetParticle(411)->Mass();
437 
438  Int_t nBinsMass = (Int_t) (fMassRange/fMassBinSize+0.5);
439  if (nBinsMass%2==1) nBinsMass++;
440  Float_t minMass = massDs - 0.5*nBinsMass*fMassBinSize;
441  Float_t maxMass = massDs + 0.5*nBinsMass*fMassBinSize;
442  Float_t minMassDplus = massDplus - 0.5*nBinsMass*fMassBinSize;
443  Float_t maxMassDplus = massDplus + 0.5*nBinsMass*fMassBinSize;
444 
445 
446  TString titleHisto(""), nameHisto("");
447  TString suffixType[5] = { "kCan", "kPid", "mcSig", "mcBackg", "mcDplusRefl" };
448  TString titleType[5] = { "kCandidate", "kCandidate + kPID", "MC Signal", "MC Background", "MC D^{+} reflection" };
449 
450 
451 
452  for (Int_t itype=0; itype<5; itype++) {
453 
454  // - Do not create MC histograms if fReadMC is not true
455  if (itype>=2 && !fReadMC) break;
456 
457 
458  nameHisto = Form("fHisInvMassDs_%s", suffixType[itype].Data());
459  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); M_{inv}(K^{0}_{S}K^{#pm}); Entries", titleType[itype].Data());
460  fHisInvMassDs[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, nBinsMass, minMass, maxMass);
461 
462  nameHisto = Form("fHisInvMassDplus_%s", suffixType[itype].Data());
463  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); M_{inv}(K^{0}_{S}#pi^{#pm}); Entries", titleType[itype].Data());
464  fHisInvMassDplus[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, nBinsMass, minMassDplus, maxMassDplus);
465 
466  nameHisto = Form("fHisInvMassK0s_%s", suffixType[itype].Data());
467  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); M_{inv}(#pi^{+},#pi^{-}); Entries", titleType[itype].Data());
468  fHisInvMassK0s[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 200, 0.4, 0.6);
469 
470  nameHisto = Form("fHisPtK0s_%s", suffixType[itype].Data());
471  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); #it{p}_{T}(K^{0}_{S}) (GeV/#it{c}); Entries", titleType[itype].Data());
472  fHisPtK0s[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 160, 0., 40.);
473 
474  nameHisto = Form("fHisPtBachelor_%s", suffixType[itype].Data());
475  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); #it{p}_{T}(K^{+}) (GeV/#it{c}); Entries", titleType[itype].Data());
476  fHisPtBachelor[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 160, 0., 40.);
477 
478  nameHisto = Form("fHisImpParK0s_%s", suffixType[itype].Data());
479  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); #it{d_{0}}(K^{0}_{S}) (cm); Entries", titleType[itype].Data());
480  fHisImpParK0s[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 1000, -5., 5.);
481 
482  nameHisto = Form("fHisImpParBach_%s", suffixType[itype].Data());
483  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); #it{d_{0}}(K^{+}) (cm); Entries", titleType[itype].Data());
484  fHisImpParBach[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 1000, -5., 5.);
485 
486  nameHisto = Form("fHisCTauK0s_%s", suffixType[itype].Data());
487  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); #it{c#times#tau}(K^{0}_{S}); Entries", titleType[itype].Data());
488  fHisCTauK0s[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 160, 0., 40.);
489 
490  nameHisto = Form("fHisCosPointingDs_%s", suffixType[itype].Data());
491  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); cos(#it{#theta_{p}})(D_{s}^{+}); Entries", titleType[itype].Data());
492  fHisCosPointingDs[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 400, 0.0, 1.);
493 
494  nameHisto = Form("fHisCosPointingXYDs_%s", suffixType[itype].Data());
495  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); cos(#it{#theta_{p}^{xy}})(D_{s}^{+}); Entries", titleType[itype].Data());
496  fHisCosPointingXYDs[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 800, -1.0, 1.);
497 
498  nameHisto = Form("fHisCosThetaStarK0s_%s", suffixType[itype].Data());
499  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); cos(#it{#theta^{*}})(K_{S}^{0}); Entries", titleType[itype].Data());
500  fHisCosThetaStarK0s[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 800, -1., 1.);
501 
502  nameHisto = Form("fHisCosThetaStarBach_%s", suffixType[itype].Data());
503  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); cos(#it{#theta^{*}})(K^{+}); Entries", titleType[itype].Data());
504  fHisCosThetaStarBach[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 800, -1., 1.);
505 
506  nameHisto = Form("fHisDCAK0sBach_%s", suffixType[itype].Data());
507  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); DCA(K^{0}_{S},K^{+}) (cm); Entries", titleType[itype].Data());
508  fHisDCAK0sBach[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 500, 0., 10.);
509 
510  nameHisto = Form("fHisDecayLxyDs_%s", suffixType[itype].Data());
511  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); #it{L_{xy}}(D^{+}_{s}) (cm); Entries", titleType[itype].Data());
512  fHisDecayLxyDs[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 400, 0., 1.0);
513 
514  nameHisto = Form("fHisNormDecayLxyDs_%s", suffixType[itype].Data());
515  titleHisto = Form("%s; #it{p}_{T}(D^{+}_{s}) (GeV/#it{c}); #it{L_{xy}}/#it{#sigma_{xy}}(D^{+}_{s}) (cm); Entries", titleType[itype].Data());
516  fHisNormDecayLxyDs[itype] = new TH2F(nameHisto.Data(), titleHisto.Data(), fNPtBins, ptBinsRange, 200, 0., 40.);
517 
518 
519  TList *theOutput = 0;
520  if (itype==0) theOutput = fOutputCand;
521  else if (itype==1) theOutput = fOutputPID;
522  else theOutput = fOutputMC;
523 
524  theOutput->Add(fHisInvMassDs[itype]);
525  theOutput->Add(fHisInvMassDplus[itype]);
526  theOutput->Add(fHisInvMassK0s[itype]);
527  theOutput->Add(fHisPtK0s[itype]);
528  theOutput->Add(fHisPtBachelor[itype]);
529  theOutput->Add(fHisImpParK0s[itype]);
530  theOutput->Add(fHisImpParBach[itype]);
531  theOutput->Add(fHisCTauK0s[itype]);
532  theOutput->Add(fHisCosPointingDs[itype]);
533  theOutput->Add(fHisCosPointingXYDs[itype]);
534  theOutput->Add(fHisCosThetaStarK0s[itype]);
535  theOutput->Add(fHisCosThetaStarBach[itype]);
536  theOutput->Add(fHisDCAK0sBach[itype]);
537  theOutput->Add(fHisDecayLxyDs[itype]);
538  theOutput->Add(fHisNormDecayLxyDs[itype]);
539 
540  } // end for
541 
542 
543 
544  } else {
545 
546  //---------------------------------------------------------------------------
547  // - Output slot #4: TNtuple for candidates on data
548  //---------------------------------------------------------------------------
549  TString listVar("K0sInvMass"); // 0 V0: invariant mass (pi+pi-)
550  listVar.Append(":K0sPt"); // 1 V0: Pt
551  listVar.Append(":K0sRxy"); // 2 V0: decay length XY
552  listVar.Append(":K0sCTau"); // 3 V0: pseudo proper decay length (L*m/p)
553  listVar.Append(":K0sCosPA"); // 4 V0: cosine pointing angle
554  listVar.Append(":K0sd0"); // 5 V0: impact parameter
555  listVar.Append(":K0sd0DauPos"); // 6 V0: impact parameter positive daughter
556  listVar.Append(":K0sd0DauNeg"); // 7 V0: impact parameter negative daughter
557  listVar.Append(":K0sDCAdau"); // 8 V0: DCA prong-to-prong (pi+,pi-)
558  listVar.Append(":K0sPtDauPos"); // 9 V0: Pt positive daughters
559  listVar.Append(":K0sPtDauNeg"); // 10 V0: Pt positive daughters
560  listVar.Append(":BachPt"); // 11 Bach: Pt
561  listVar.Append(":Bachd0"); // 12 Bach: impact parameter
562  listVar.Append(":BachNsigmaTPCpion"); // 13 Bach: N sigmaTPC pion
563  listVar.Append(":BachNsigmaTPCkaon"); // 14 Bach: N sigmaTPC kaon
564  listVar.Append(":BachNsigmaTPCproton"); // 15 Bach: N sigmaTPC proton
565  listVar.Append(":BachNsigmaTOFpion"); // 16 Bach: N sigmaTOF pion
566  listVar.Append(":BachNsigmaTOFkaon"); // 17 Bach: N sigmaTOF kaon
567  listVar.Append(":BachNsigmaTOFproton"); // 18 Bach: N sigmaTOF proton
568  listVar.Append(":CanInvMassDs"); // 19 Ds: invariant mass (K0sK+)
569  listVar.Append(":CanInvMassDplus"); // 20 Ds: invariant mass (K0spi+)
570  listVar.Append(":CanPt"); // 21 Ds: Pt
571  listVar.Append(":CanDCAProngToProng"); // 22 Ds: DCA prong-to-prong (K0s,K)
572  listVar.Append(":CanCosThetaStarK0s"); // 23 Ds: cosine theta* K0s
573  listVar.Append(":CanCosThetaStarBach"); // 24 Ds: cosine theta* K+
574  listVar.Append(":CanCosPA"); // 25 Ds: cosine pointing angle
575  listVar.Append(":CanCosPAxy"); // 26 Ds: cosine pointing angle XY
576  listVar.Append(":CanDLengthXY"); // 27 Ds: decay length XY
577  listVar.Append(":CanNormDLengthXY"); // 28 Ds: normalised decay length XY
578  listVar.Append(":CanDLength3D"); // 29 Ds: decay length 3D
579  listVar.Append(":CanNormDLength3D"); // 30 Ds: normalised decay length 3D
580  listVar.Append(":CanSigmaVtx"); // 31 Ds: sigma vertex
581  listVar.Append(":CanNormTopoBach"); // 32 Ds: difference between measured and expected normalised d0(K+)
582  listVar.Append(":CanNormTopoK0s"); // 33 Ds: difference between measured and expected normalised d0(K0s)
583  if (fReadMC) {
584  listVar.Append(":mcTruth"); // 34 MC: MC informations: 10 (signal), 20 (D+ reflection),
585  // -1 (backg w/ true K0S), -2 (backg w/ false K0S)
586  }
587  // listVar.Append(":CanCosThetaRFrame"); // 32 Ds: cosine of angle between (K0s-bach) in the Ds rest frame
588 
589 
590  fOutputNtuple = new TNtuple(GetOutputSlot(4)->GetContainer()->GetName(), "TNtuple for reconstructed candidates on real data", listVar.Data());
591 
592  }
593 
594 
595 
596  PostData(2, fCounter);
597  PostData(3, fOutputSele);
598  if (!fFillNtuple) {
599  PostData(4, fOutputCand);
600  PostData(5, fOutputPID);
601  if (fReadMC)
602  PostData(6, fOutputMC);
603  } else {
604  PostData(4, fOutputNtuple);
605  }
606 
607  return;
608 }
609 
610 
611 
612 
613 
614 //__________________________________________________________________________
616 {
620  AliDebug(1, "AliAnalysisTaskSEDstoK0sK::UserExec()");
621 
622 
623  //---------------------------------------------------------------------------
624  // - Load the event and the CascadesHF branch
625  //---------------------------------------------------------------------------
626 
627  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
628  fHisNEvents->Fill(0); // all events
629 
630 
631  if(fAODProtection>=0){
632  // Protection against different number of events in the AOD and deltaAOD
633  // In case of discrepancy the event is rejected.
634  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
635  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
636  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
637  fHisNEvents->Fill(2);
638  return;
639  }
640  fHisNEvents->Fill(1);
641  }
642 
643 
644  TClonesArray *arrayCascadesHF = 0;
645  if(!aod && AODEvent() && IsStandardAOD()) {
646  // In case there is an AOD handler writing a standard AOD, use the AOD
647  // event in memory rather than the input (ESD) event.
648  aod = dynamic_cast<AliAODEvent*> (AODEvent());
649  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
650  // have to taken from the AOD event hold by the AliAODExtension
651  AliAODHandler *aodHandler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
652  if (aodHandler->GetExtensions()) {
653  AliAODExtension *ext = (AliAODExtension*) aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
654  AliAODEvent *aodFromExt = ext->GetAOD();
655  arrayCascadesHF = (TClonesArray*) aodFromExt->GetList()->FindObject("CascadesHF");
656  }
657  } else if (aod) {
658  arrayCascadesHF = (TClonesArray*) aod->GetList()->FindObject("CascadesHF");
659  }
660 
661  if (!aod || !arrayCascadesHF) {
662  AliWarning("AliAnalysisTaskSEDstoK0sK::UserExec: CascadesHF branch not found!");
663  return;
664  }
665 
666 
667  // Load MC particles
668  TClonesArray *mcArray = 0;
669  AliAODMCHeader *mcHeader = 0;
670 
671  if (fReadMC) {
672  mcArray = (TClonesArray*) aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
673  if (!mcArray) {
674  AliWarning("AliAnalysisTaskSEDstoK0sK::UserExec: MC particles branch not found!");
675  return;
676  }
677  mcHeader = (AliAODMCHeader*) aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
678  if (!mcHeader) {
679  AliWarning("AliAnalysisTaskSEDstoK0sK::UserExec: MC header branch not found!");
680  return;
681  }
682 
683  // FillMCAcceptanceHistos has to take place here
684  }
685 
686 
687 
688  //---------------------------------------------------------------------------
689  // - Event selection
690  //---------------------------------------------------------------------------
691 
692  // Fix for temporary bug in ESDfilter
693  // The AODs with null vertex pointer didn't pass the PhysSel
694  if (!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
695 
696  fHisNEvents->Fill(3);
698 
699 
700  Float_t nTracks = aod->GetNumberOfTracks();
702 
703  fHisCentrality[0]->Fill(evCentr);
704  fHisCentralityVsMult[0]->Fill(nTracks, evCentr);
705 
706 
707  if ( !(fAnalysisCuts->IsEventSelected(aod)) ) {
708 
709  if (fAnalysisCuts->GetWhyRejection() == 5) fHisNEvents->Fill(4); // Rejection: Trigger
710  if (fAnalysisCuts->GetWhyRejection() == 0) fHisNEvents->Fill(5); // Rejection: Vertex reco
711  if (fAnalysisCuts->GetWhyRejection() == 1) fHisNEvents->Fill(6); // Rejection: Pileup
712  if (fAnalysisCuts->GetWhyRejection() == 2) fHisNEvents->Fill(7); // Rejection: Centrality
713  if (fAnalysisCuts->GetWhyRejection() == 6) fHisNEvents->Fill(8); // Rejection: VtxZ fiducial acc.
714  if (fAnalysisCuts->GetWhyRejection() == 7) fHisNEvents->Fill(9); // Rejection: Physics selection
715 
716  if (fAnalysisCuts->GetWhyRejection() == 2) {
717  fHisCentrality[2]->Fill(evCentr);
718  fHisCentralityVsMult[2]->Fill(nTracks, evCentr);
719  }
720 
721  return;
722  }
723 
724 
725  fHisNEvents->Fill(10);
726  fHisCentrality[1]->Fill(evCentr);
727  fHisCentralityVsMult[1]->Fill(nTracks, evCentr);
728 
729 
730  Int_t nCascades = arrayCascadesHF->GetEntriesFast();
731  AliDebug(1, Form("AliAnalysisTaskSEDstoK0sK::UserExec: %d Cascades found for this event!", nCascades));
732 
733 
734 
735 
736  //---------------------------------------------------------------------------
737  // - Analysis: Loop over the candidates
738  //---------------------------------------------------------------------------
739 
740  // vHF object is needed to call the method that refills the missing info of the candidates
741  // if they have been deleted in dAOD reconstruction phase
742  // in order to reduce the size of the file
744 
745 
746 
747  for (Int_t iCasc=0; iCasc<nCascades; iCasc++)
748  {
749 
750  //-----------------------------------------
751  // - Get candidates and fill missing info
752  //-----------------------------------------
753  AliAODRecoCascadeHF *dCan = (AliAODRecoCascadeHF*) arrayCascadesHF->UncheckedAt(iCasc);
754 
755  if (fUseSelectionBit) {
756  if ( !(dCan->CheckCascadeFlags(AliRDHFCuts::kDstoK0sCuts)) ) continue;
757  }
758 
759  fHisNEvents->Fill(11);
760 
761 
762  // Fill the data members of the candidate only if they are empty.
763  if (!(vHF->FillRecoCasc(aod, dCan, kFALSE, kTRUE))) {
764  fHisNEvents->Fill(12);
765  continue;
766  }
767 
768 
769  // Reconstruct secondary vertex if it does not already exist (standard dAOD)
770  if (dCan->GetIsFilled()==1) {
771  if (!(vHF->RecoSecondaryVertexForCascades(aod, dCan))) {
772  fHisNEvents->Fill(13);
773  continue;
774  }
775  }
776 
777  fHisNEvents->Fill(14);
778 
779 
780 
781  //-----------------------------------------
782  // - Select candidates step by step
783  //-----------------------------------------
784  Double_t ptCand = dCan->Pt();
785  Double_t rapidity = dCan->Y(431);
786 
787  fHisRapidity->Fill(ptCand, rapidity);
788 
789 
790  if (!(fAnalysisCuts->IsInFiducialAcceptance(ptCand, rapidity))) continue;
791  fHisNEvents->Fill(15);
792 
793 
794  if (!(fAnalysisCuts->IsSelected(dCan, AliRDHFCuts::kTracks, aod))) continue;
795  fHisNEvents->Fill(16);
796 
797  // FillTheTree(dCan, mcArray);
798 
799 
800  if (!(fAnalysisCuts->IsSelected(dCan, AliRDHFCuts::kCandidate, aod))) continue;
801  fHisNEvents->Fill(17);
802 
805 
806 
807  if (!(fAnalysisCuts->IsSelected(dCan, AliRDHFCuts::kPID, aod))) continue;
808  fHisNEvents->Fill(18);
809 
810  fHisRapiditySel->Fill(ptCand, rapidity);
811 
812  FillHistogramsVar(dCan, AliRDHFCuts::kPID, mcArray);
813  FillHistogramsPID(dCan, AliRDHFCuts::kPID, mcArray);
814 
815 
816  FillTheTree(dCan, aod, mcArray);
817 
818 
819  } // end loop cascades
820 
821 
822 
823  delete vHF;
824 
825 
826  PostData(2, fCounter);
827  PostData(3, fOutputSele);
828  if (!fFillNtuple) {
829  PostData(4, fOutputCand);
830  PostData(5, fOutputPID);
831  if (fReadMC)
832  PostData(6, fOutputMC);
833  } else {
834  PostData(4, fOutputNtuple);
835  }
836 
837  return;
838 }
839 
840 
841 
842 
843 //__________________________________________________________________________
845 {
849  AliDebug(1, "AliAnalysisTaskSEDstoK0sK::Terminate()");
850 
851 
852  // Output slot #3
853  fOutputSele = dynamic_cast<TList*> (GetOutputData(3));
854  if (!fOutputSele) {
855  AliWarning("ERROR: fOutputSele not available");
856  return;
857  }
858 
859 
860  if (!fFillNtuple) {
861 
862  // Output slot #4
863  fOutputCand = dynamic_cast<TList*> (GetOutputData(4));
864  if (!fOutputCand) {
865  AliWarning("ERROR: fOutputCand not available");
866  return;
867  }
868 
869  // Output slot #5
870  fOutputPID = dynamic_cast<TList*> (GetOutputData(5));
871  if (!fOutputPID) {
872  AliWarning("ERROR: fOutputPID not available");
873  return;
874  }
875 
876  if (fReadMC) {
877  // Output slot #6
878  fOutputMC = dynamic_cast<TList*> (GetOutputData(6));
879  if (!fOutputMC) {
880  AliWarning("ERROR: fOutputMC not available");
881  return;
882  }
883  }
884 
885  } else {
886 
887  // Output slot #4
888  fOutputNtuple = dynamic_cast<TNtuple*> (GetOutputData(4));
889  if (!fOutputNtuple) {
890  AliWarning("ERROR: fOutputNtuple not available");
891  return;
892  }
893 
894  }
895 
896  return;
897 }
898 
899 
900 
901 
902 //__________________________________________________________________________
904 {
908 
909 
910  if (nBins<=kMaxPtBins) {
911 
912  fNPtBins = nBins;
913  for (Int_t ipt=0; ipt<fNPtBins+1; ipt++) fPtLimits[ipt] = limitsPt[ipt];
914  for (Int_t ipt=fNPtBins+1; ipt<kMaxPtBins+1; ipt++) fPtLimits[ipt] = 99999999.;
915 
916  } else {
917 
918  // More Pt bins than allowed by the class header
919  AliWarning(Form("Maximum number of Pt bins reached: %d", kMaxPtBins));
921  fPtLimits[0] = 2.;
922  fPtLimits[1] = 4.;
923  fPtLimits[2] = 6.;
924  fPtLimits[3] = 8.;
925  for (Int_t ipt=4; ipt<kMaxPtBins+1; ipt++) fPtLimits[ipt] = 99999999.;
926 
927  }
928 
929  if (fDebug>1) {
930  AliInfo(Form("Number of Pt bins for the analysis: %d", fNPtBins));
931  for (Int_t ipt=0; ipt<fNPtBins; ipt++) {
932  AliInfo(Form(" Pt bin %02d: [%8.1f, %8.1f]", ipt, fPtLimits[ipt], fPtLimits[ipt+1]));
933  }
934  }
935 
936  return;
937 }
938 
939 
940 
941 
942 //__________________________________________________________________________
944 {
948 
949 
950  if (nCuts==kNTupleVars) {
951 
952  for (Int_t ivar=0; ivar<kNTupleVars; ivar++) {
953  fCutsMinTupleVars[ivar] = minCuts[ivar];
954  fCutsMaxTupleVars[ivar] = maxCuts[ivar];
955  }
956 
957  } else {
958 
959  // Different number of cuts than the number of variables in the tuple
960  AliWarning(Form("%d cuts are set for %d variables in the tuple: default initialisation", nCuts, kNTupleVars));
961 
962  for (Int_t ivar=0; ivar<kNTupleVars; ivar++) {
963  fCutsMinTupleVars[ivar] = -999.;
964  fCutsMaxTupleVars[ivar] = 999.;
965  }
966 
967  }
968 
969 
970  return;
971 }
972 
973 
974 
975 
976 //__________________________________________________________________________
978 {
989  AliDebug(1, "AliAnalysisTaskSEDstoK0sK::FillHistogramsVar()");
990 
991 
992  if (fFillNtuple) return;
993 
994 
995  //---------------------------------------------------------------------------
996  // - Monte Carlo: check the candidate at MC level
997  //---------------------------------------------------------------------------
998  Bool_t okMCsignal = kFALSE;
999  Bool_t okMCbackg = kFALSE;
1000  Bool_t okMCreflec = kFALSE;
1001 
1002  if (selFlag==AliRDHFCuts::kPID && fReadMC) {
1003  if (!mcArray) return;
1004 
1005  if (MatchToMCDstoK0sKSignal(dCan, mcArray)>=0) {
1006  AliDebug(2, "This candidate matches a MC signal of Ds+ -> K0S + K+");
1007  okMCsignal = kTRUE;
1008  } else if (MatchToMCDplustoK0spiSignal(dCan, mcArray)>=0) {
1009  AliDebug(2, "This candidate matches a MC signal of D+ -> K0S + pi+");
1010  okMCreflec = kTRUE;
1011  } else {
1012  AliDebug(2, "This candidate is a MC background");
1013  okMCbackg = kTRUE;
1014  }
1015  }
1016 
1017 
1018 
1019  //---------------------------------------------------------------------------
1020  // - Fill the histograms depending on the selection level
1021  // - if 'kCandidate' level, fill only index=0 (histograms 'kCand')
1022  // - if 'kPID' level and not reading MC, fill only index=1 (histograms 'kPID')
1023  // - if 'kPID' level and reading MC, fill index=1 ('kPID'), 2 ('mcSig'), 3 ('mcBackg'), 4 ('mcDplusRefl')
1024  //---------------------------------------------------------------------------
1025  Float_t ptCand = dCan->Pt();
1026 
1027 
1028  for (Int_t index=0; index<5; index++) {
1029 
1030  if (index==0) {
1031  if (selFlag != AliRDHFCuts::kCandidate) continue;
1032  } else if (index==1) {
1033  if (selFlag != AliRDHFCuts::kPID) continue;
1034  } else if (index==2) {
1035  if (!okMCsignal) continue;
1036  } else if (index==3) {
1037  if (!okMCbackg) continue;
1038  } else if (index==4) {
1039  if (!okMCreflec) continue;
1040  }
1041 
1042  fHisInvMassDs[index] ->Fill( ptCand, dCan->InvMassDstoK0sK() );
1043  fHisInvMassDplus[index] ->Fill( ptCand, dCan->InvMassDplustoK0spi() );
1044  fHisInvMassK0s[index] ->Fill( ptCand, dynamic_cast<AliAODv0*>(dCan->Getv0())->MassK0Short());
1045  fHisPtK0s[index] ->Fill( ptCand, dCan->PtProng(1) );
1046  fHisPtBachelor[index] ->Fill( ptCand, dCan->PtProng(0) );
1047  fHisImpParK0s[index] ->Fill( ptCand, dCan->Getd0Prong(1) );
1048  fHisImpParBach[index] ->Fill( ptCand, dCan->Getd0Prong(0) );
1049  fHisCTauK0s[index] ->Fill( ptCand, dynamic_cast<AliAODv0*>(dCan->Getv0())->Ct(310, dCan->GetOwnPrimaryVtx()));
1050  fHisCosPointingDs[index] ->Fill( ptCand, dCan->CosPointingAngle() );
1051  fHisCosPointingXYDs[index] ->Fill( ptCand, dCan->CosPointingAngleXY() );
1052  fHisCosThetaStarK0s[index] ->Fill( ptCand, dCan->CosThetaStar(1, 431, 321, 310) );
1053  fHisCosThetaStarBach[index] ->Fill( ptCand, dCan->CosThetaStar(0, 431, 321, 310) );
1054  fHisDCAK0sBach[index] ->Fill( ptCand, dCan->GetDCA() );
1055  fHisDecayLxyDs[index] ->Fill( ptCand, dCan->DecayLengthXY() );
1056  fHisNormDecayLxyDs[index] ->Fill( ptCand, dCan->NormalizedDecayLengthXY() );
1057 
1058  }
1059 
1060 
1061  return;
1062 }
1063 
1064 
1065 
1066 
1067 //__________________________________________________________________________
1069 {
1073  AliDebug(2, "AliAnalysisTaskSEDstoK0sK::FillHistogramsPID()");
1074 
1075 
1076  Double_t nSigmaTPC = -999;
1077  Double_t nSigmaTOF = -999.;
1078  Float_t ptCand = dCan->Pt();
1079 
1080  if (fAnalysisCuts->GetIsUsePID()) {
1082  if (pidhf->CheckTPCPIDStatus(dCan->GetBachelor())) pidhf->GetnSigmaTPC(dCan->GetBachelor(), AliPID::kKaon, nSigmaTPC);
1083  if (pidhf->CheckTOFPIDStatus(dCan->GetBachelor())) pidhf->GetnSigmaTOF(dCan->GetBachelor(), AliPID::kKaon, nSigmaTOF);
1084  }
1085 
1086 
1087  if (selFlag==AliRDHFCuts::kCandidate) {
1088 
1089  //---------------------------------------------------------------------------
1090  // - Topological selection levels
1091  //---------------------------------------------------------------------------
1092  fHisPidTPCKaonVsPt->Fill(ptCand, nSigmaTPC);
1093  fHisPidTOFKaonVsPt->Fill(ptCand, nSigmaTOF);
1094  fHisPidTPCTOFKaon->Fill(nSigmaTPC, nSigmaTOF);
1095 
1096  } else if (selFlag==AliRDHFCuts::kPID || selFlag==AliRDHFCuts::kAll) {
1097 
1098  //---------------------------------------------------------------------------
1099  // - Topological+PID selection levels
1100  //---------------------------------------------------------------------------
1101  fHisPidTPCKaonVsPtSel->Fill(ptCand, nSigmaTPC);
1102  fHisPidTOFKaonVsPtSel->Fill(ptCand, nSigmaTOF);
1103  fHisPidTPCTOFKaonSel->Fill(nSigmaTPC, nSigmaTOF);
1104 
1105  //-----------------------------------------
1106  // - Check the candidate at MC level
1107  //-----------------------------------------
1108  if (fReadMC && mcArray) {
1109  if (MatchToMCDstoK0sKSignal(dCan, mcArray)>=0) fHisNEvents->Fill(19);
1110  else if (MatchToMCDplustoK0spiSignal(dCan, mcArray)>=0) fHisNEvents->Fill(20);
1111  }
1112 
1113  }
1114 
1115 
1116  return;
1117 }
1118 
1119 
1120 
1121 
1122 //__________________________________________________________________________
1124 {
1128  AliDebug(2, "AliAnalysisTaskSEDstoK0sK::FillTheTree()");
1129 
1130 
1131  if (!fFillNtuple) return;
1132 
1133  AliAODv0 *v0 = (AliAODv0*) dCan->Getv0();
1134  if (!v0) return;
1135 
1136 
1137  //---------------------------------------------------------------------------
1138  // - Get PID informations if possible
1139  //---------------------------------------------------------------------------
1140  Double_t nSigmaTPCpion = -999;
1141  Double_t nSigmaTPCkaon = -999;
1142  Double_t nSigmaTPCproton = -999;
1143  Double_t nSigmaTOFpion = -999;
1144  Double_t nSigmaTOFkaon = -999;
1145  Double_t nSigmaTOFproton = -999;
1146  if (fAnalysisCuts->GetIsUsePID()) {
1148  fAnalysisCuts->GetPidHF()->GetnSigmaTPC(dCan->GetBachelor(), AliPID::kPion, nSigmaTPCpion);
1149  fAnalysisCuts->GetPidHF()->GetnSigmaTPC(dCan->GetBachelor(), AliPID::kKaon, nSigmaTPCkaon);
1150  fAnalysisCuts->GetPidHF()->GetnSigmaTPC(dCan->GetBachelor(), AliPID::kProton, nSigmaTPCproton);
1151  }
1153  fAnalysisCuts->GetPidHF()->GetnSigmaTOF(dCan->GetBachelor(), AliPID::kPion, nSigmaTOFpion);
1154  fAnalysisCuts->GetPidHF()->GetnSigmaTOF(dCan->GetBachelor(), AliPID::kKaon, nSigmaTOFkaon);
1155  fAnalysisCuts->GetPidHF()->GetnSigmaTOF(dCan->GetBachelor(), AliPID::kProton, nSigmaTOFproton);
1156  }
1157  }
1158 
1159 
1160  //---------------------------------------------------------------------------
1161  // - Apply cuts to the tuple variables
1162  //---------------------------------------------------------------------------
1163  Double_t diffIP[2] = {0.}, errdiffIP[2] = {1.};
1164  dCan->Getd0MeasMinusExpProng(0, aod->GetMagneticField(), diffIP[0], errdiffIP[0]);
1165  dCan->Getd0MeasMinusExpProng(1, aod->GetMagneticField(), diffIP[1], errdiffIP[1]);
1166 
1167 
1168  const Int_t nVarTuple = fReadMC ? kNTupleVarsMC : kNTupleVars;
1169 
1170  Float_t variableNtuple[nVarTuple];
1171  variableNtuple[ 0] = v0->MassK0Short(); // "K0sInvMass"
1172  variableNtuple[ 1] = dCan->PtProng(1); // "K0sPt"
1173  variableNtuple[ 2] = v0->RadiusV0(); // "K0sRxy"
1174  variableNtuple[ 3] = v0->Ct(310, dCan->GetOwnPrimaryVtx()); // "K0sCTau"
1175  variableNtuple[ 4] = dCan->CosV0PointingAngle(); // "K0sCosPA"
1176  variableNtuple[ 5] = dCan->Getd0Prong(1); // "K0sd0"
1177  variableNtuple[ 6] = v0->DcaPosToPrimVertex(); // "K0sd0DauPos"
1178  variableNtuple[ 7] = v0->DcaNegToPrimVertex(); // "K0sd0DauNeg"
1179  variableNtuple[ 8] = v0->DcaV0Daughters(); // "K0sDCAdau"
1180  variableNtuple[ 9] = v0->PtProng(0); // "K0sPtDauPos"
1181  variableNtuple[10] = v0->PtProng(1); // "K0sPtDauNeg"
1182  variableNtuple[11] = dCan->PtProng(0); // "BachPt"
1183  variableNtuple[12] = dCan->Getd0Prong(0); // "Bachd0"
1184  variableNtuple[13] = nSigmaTPCpion; // "BachNsigmaTPCpion"
1185  variableNtuple[14] = nSigmaTPCkaon; // "BachNsigmaTPCkaon"
1186  variableNtuple[15] = nSigmaTPCproton; // "BachNsigmaTPCproton"
1187  variableNtuple[16] = nSigmaTOFpion; // "BachNsigmaTOFpion"
1188  variableNtuple[17] = nSigmaTOFkaon; // "BachNsigmaTOFkaon"
1189  variableNtuple[18] = nSigmaTOFproton; // "BachNsigmaTOFproton"
1190  variableNtuple[19] = dCan->InvMassDstoK0sK(); // "CanInvMassDs"
1191  variableNtuple[20] = dCan->InvMassDplustoK0spi(); // "CanInvMassDplus"
1192  variableNtuple[21] = dCan->Pt(); // "CanPt"
1193  variableNtuple[22] = dCan->GetDCA(); // "CanDCAProngToProng"
1194  variableNtuple[23] = dCan->CosThetaStar(1, 431, 321, 310); // "CanCosThetaStarK0s"
1195  variableNtuple[24] = dCan->CosThetaStar(0, 431, 321, 310); // "CanCosThetaStarBach"
1196  variableNtuple[25] = dCan->CosPointingAngle(); // "CanCosPA"
1197  variableNtuple[26] = dCan->CosPointingAngleXY(); // "CanCosPAxy"
1198  variableNtuple[27] = dCan->DecayLengthXY(); // "CanDLengthXY"
1199  variableNtuple[28] = dCan->NormalizedDecayLengthXY(); // "CanNormDLengthXY"
1200  variableNtuple[29] = dCan->DecayLength(); // "CanDLength3D"
1201  variableNtuple[30] = dCan->NormalizedDecayLength(); // "CanNormDLength3D"
1202  variableNtuple[31] = ComputeSigmaVert(aod, dCan); // "CanSigmaVtx"
1203  variableNtuple[32] = diffIP[0]/errdiffIP[0]; // "CanNormTopoBach"
1204  variableNtuple[33] = diffIP[1]/errdiffIP[1]; // "CanNormTopoK0s"
1205  // variableNtuple[32] = CosThetaK0sBachRFrame(dCan); // "CanCosThetaRFrame"
1206 
1207 
1208  for (Int_t ivar=0; ivar<kNTupleVars; ivar++) {
1209  if ( (variableNtuple[ivar]<fCutsMinTupleVars[ivar]) || (variableNtuple[ivar]>fCutsMaxTupleVars[ivar]) ) {
1210  AliDebug(2, "One of the tuple variable does not pass the selection cut");
1211  return;
1212  }
1213  }
1214 
1215 
1216  //---------------------------------------------------------------------------
1217  // - Monte Carlo: check the candidate at MC level
1218  //---------------------------------------------------------------------------
1219  if (fReadMC) {
1220  Bool_t okMCsignal = (MatchToMCDstoK0sKSignal(dCan, mcArray)>=0) ? kTRUE : kFALSE;
1221  Bool_t okMCreflec = (MatchToMCDplustoK0spiSignal(dCan, mcArray)>=0) ? kTRUE : kFALSE;
1222 
1223  variableNtuple[34] = okMCsignal ? 10 // "CanNormDLengthXY"
1224  : okMCreflec ? 20
1225  : -2;
1226 
1227  if (!okMCsignal && !okMCreflec) {
1228  AliDebug(2, "The candidate is a combinatorial background, check if the V0 is a true K0S");
1229  Int_t pdgDgK0stoPions[2] = {211, 211};
1230  if (v0->MatchToMC(310, mcArray, 2, pdgDgK0stoPions)>=0) variableNtuple[29]++;
1231  }
1232 
1233  } // end if fReadMC
1234 
1235 
1236  //---------------------------------------------------------------------------
1237  // - Fill the tuple
1238  //---------------------------------------------------------------------------
1239  fOutputNtuple->Fill(variableNtuple);
1240 
1241 
1242  return;
1243 }
1244 
1245 
1246 
1247 
1248 //__________________________________________________________________________
1250 {
1255 
1256 
1257  // - Get bachelor and V0 from the cascade
1258  AliESDtrack *esdB = new AliESDtrack(dCan->GetBachelor());
1259 
1260  AliNeutralTrackParam *trackV0 = 0;
1261 
1262  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(dCan->Getv0());
1263  if (trackVV0) {
1264  trackV0 = new AliNeutralTrackParam(trackVV0);
1265  }
1266 
1267 
1268  if (!esdB || !trackV0) {
1269  if (esdB) delete esdB;
1270  if (trackV0) delete trackV0;
1271  return -999.;
1272  }
1273 
1274 
1275  // - Build ESD primary vertex
1276  AliVertexerTracks vertexer(aod->GetMagneticField());
1277 
1278  Double_t pos[3], cov[6];
1279  AliAODVertex *aodV = (AliAODVertex*) aod->GetPrimaryVertex();
1280  aodV->GetXYZ(pos);
1281  aodV->GetCovarianceMatrix(cov);
1282  Double_t chi2 = aodV->GetChi2();
1283  Int_t nContr = aodV->GetNContributors();
1284 
1285  AliESDVertex vprim(pos, cov, chi2, nContr);
1286  vertexer.SetVtxStart(&vprim);
1287 
1288 
1289  // - Assign TObjArray of daughters
1290  TObjArray twoTrackArray(2);
1291  twoTrackArray.AddAt(esdB, 0);
1292  twoTrackArray.AddAt(trackV0, 1);
1293 
1294 
1295  // - Build ESD secondary vertex and get dispersion
1296  AliESDVertex *secVert = vertexer.VertexForSelectedESDTracks(&twoTrackArray, kFALSE, kTRUE, kFALSE);
1297  Double_t disp = secVert->GetDispersion();
1298 
1299 
1300  twoTrackArray.Delete();
1301  delete secVert;
1302 
1303  return disp;
1304 }
1305 
1306 
1307 
1308 
1309 //__________________________________________________________________________
1311 {
1316 
1317 
1318  // - Get bachelor and V0 from the cascade
1319  AliAODv0 *v0 = dynamic_cast<AliAODv0*> (dCan->Getv0());
1320  AliAODTrack *bach = dynamic_cast<AliAODTrack*> (dCan->GetBachelor());
1321  if (!v0 || !bach) {
1322  return -999;
1323  }
1324 
1325 
1326  // - Get cascade informations necessary to compute the rest frame
1327  Double_t eCasc = dCan->EProng(0, 321) + dCan->EProng(1, 310);
1328  Double_t pxCasc = dCan->PxProng(0) + dCan->PxProng(1);
1329  Double_t pyCasc = dCan->PyProng(0) + dCan->PyProng(1);
1330  Double_t pzCasc = dCan->PzProng(0) + dCan->PzProng(1);
1331  Double_t bxCasc = pxCasc/eCasc;
1332  Double_t byCasc = pyCasc/eCasc;
1333  Double_t bzCasc = pzCasc/eCasc;
1334 
1335 
1336  // - Boost the K+ and K0s in the cascade rest frame
1337  TVector3 vecK0sCascFrame;
1338  TLorentzVector vecK0s(dCan->PxProng(1), dCan->PyProng(1), dCan->PzProng(1), dCan->EProng(1, 310));
1339  vecK0s.Boost(-bxCasc, -byCasc, -bzCasc);
1340  vecK0s.Boost(vecK0sCascFrame);
1341  vecK0sCascFrame = vecK0s.BoostVector();
1342 
1343  TVector3 vecBachCascFrame;
1344  TLorentzVector vecBach(dCan->PxProng(0), dCan->PyProng(0), dCan->PzProng(0), dCan->EProng(0, 321));
1345  vecBach.Boost(-bxCasc, -byCasc, -bzCasc);
1346  vecBach.Boost(vecBachCascFrame);
1347  vecBachCascFrame = vecBach.BoostVector();
1348 
1349 
1350  // - Compute the cosine of the angle between K+ and K0s
1351  Double_t scalProdBachK0s = vecBachCascFrame.Dot(vecK0sCascFrame);
1352  Double_t normBach = vecBachCascFrame.Mag();
1353  Double_t normK0s = vecK0sCascFrame.Mag();
1354  Double_t cosAngleCascFrame = scalProdBachK0s/(normBach*normK0s);
1355 
1356 
1357  return cosAngleCascFrame;
1358 }
1359 
Bool_t fReadMC
Bin size of invariant mass histograms (GeV)
Double_t NormalizedDecayLengthXY() const
TList * fOutputCand
! Candidate level histograms: TList sent to output slot #4
Double_t NormalizedDecayLength() const
Bool_t fUseSelectionBit
Flag for accessing MC.
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
double Double_t
Definition: External.C:58
AliRDHFCutsDstoK0sK * fAnalysisCuts
TH2F * fHisDecayLxyDs[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
Definition: External.C:236
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
Double_t InvMassDstoK0sK() const
Int_t GetnSigmaTOF(AliAODTrack *track, Int_t species, Double_t &sigma) const
Int_t fAODProtection
Maximum cut values for tuple variables.
Double_t Ct(UInt_t pdg) const
Float_t fMassRange
Limits of Pt bins.
Int_t GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &sigma) const
Bool_t CheckTOFPIDStatus(AliAODTrack *track) const
AliAODv0 * Getv0() const
Bool_t CheckTPCPIDStatus(AliAODTrack *track) const
Float_t fMassBinSize
Size range of invariant mass histograms.
static Int_t CheckMatchingAODdeltaAODevents()
Bool_t fFillNtuple
Flag for using selection bit (to select Ds flags)
Float_t ComputeSigmaVert(const AliAODEvent *aod, AliAODRecoCascadeHF *dCan) const
AliNormalizationCounter * fCounter
Cut object for Analysis on output slot #1.
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:304
Double_t CosPointingAngleXY() const
TH2F * fHisCosPointingXYDs[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
TH2F * fHisInvMassK0s[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
TList * fOutputMC
! PID level histograms: TList sent to output slot #5
TNtuple * fOutputNtuple
! TNtuple for candidates on data sent to output slot #4
Float_t fCutsMinTupleVars[kNTupleVars]
Flag for using THnSparse.
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:237
void FillHistogramsPID(AliAODRecoCascadeHF *dCan, AliRDHFCuts::ESelLevel selFlag, TClonesArray *mcArray=0)
Float_t fCutsMaxTupleVars[kNTupleVars]
Minimum cut values for tuple variables.
virtual void Terminate(Option_t *)
int Int_t
Definition: External.C:63
TH2F * fHisCosThetaStarBach[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
Int_t GetIsFilled() const
float Float_t
Definition: External.C:68
Class for cuts on AOD reconstructed Ds->K0S+K.
void SetCutsTupleVariables(Int_t nCuts, Float_t *minCuts, Float_t *maxCuts)
AliAnalysisTaskSE to produce Ds->K0S+K invariant mass spectra and THnSparse for cut optimisations...
AliAODTrack * GetBachelor() const
void FillTheTree(AliAODRecoCascadeHF *dCan, AliAODEvent *aod, TClonesArray *mcArray)
Int_t MatchToMCDplustoK0spiSignal(AliAODRecoCascadeHF *dCan, TClonesArray *mcArray)
AliAODVertex * GetOwnPrimaryVtx() const
rapidity
Definition: HFPtSpectrum.C:47
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:248
void SetPtBins(Int_t nBins, Float_t *limitsPt)
TH2F * fHisInvMassDplus[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
Double_t CosV0PointingAngle() const
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)
Float_t CosThetaK0sBachRFrame(AliAODRecoCascadeHF *dCan) const
Double_t DecayLengthXY() const
ClassImp(AliAnalysisTaskDeltaPt) AliAnalysisTaskDeltaPt
TH2F * fHisCosThetaStarK0s[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Bool_t IsEventSelected(AliVEvent *event)
TH2F * fHisPtK0s[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
TH2F * fHisRapiditySel
! Rapidity selected (kCandidate) candidates
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
TH1F * fHisCentrality[3]
! Centrality: all, selected and rejected
TList * fOutputSele
! Various histograms of selected events: TList sent to output slot #3
Double_t minMass
Float_t * GetPtBinLimits() const
Definition: AliRDHFCuts.h:238
Bool_t CheckCascadeFlags(AliRDHFCuts::ESele selFlag=AliRDHFCuts::kLctoV0Cuts)
TH2F * fHisDCAK0sBach[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
TH1F * fHisNEvents
! Counter: events and candidates
TH2F * fHisImpParBach[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
Float_t fPtLimits[kMaxPtBins+1]
Number of Pt bins.
TList * fOutputPID
! PID level histograms: TList sent to output slot #5
Bool_t GetIsUsePID() const
Definition: AliRDHFCuts.h:260
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:239
void FillHistogramsVar(AliAODRecoCascadeHF *dCan, AliRDHFCuts::ESelLevel selFlag, TClonesArray *mcArray=0)
Double_t maxMass
TH2F * fHisInvMassDs[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Bool_t RecoSecondaryVertexForCascades(AliVEvent *event, AliAODRecoCascadeHF *rc)
TH2F * fHisImpParK0s[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
TH2F * fHisCTauK0s[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
TH2F * fHisCosPointingDs[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:270
TH2F * fHisPtBachelor[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
Double_t DecayLength() const
Double_t InvMassDplustoK0spi() const
TH2F * fHisNormDecayLxyDs[5]
! (kCandidate, kPID, mcSignal, mcBackground, mcReflection)
TH2F * fHisCentralityVsMult[3]
! Centrality VS Multiplicity: all, selected and rejected
Int_t MatchToMCDstoK0sKSignal(AliAODRecoCascadeHF *dCan, TClonesArray *mcArray)
TH2F * fHisRapidity
! Rapidity all candidates