AliPhysics  b76e98e (b76e98e)
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)
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(":CanInvMassDs"); // 13 Ds: invariant mass (K0sK+)
563  listVar.Append(":CanInvMassDplus"); // 14 Ds: invariant mass (K0spi+)
564  listVar.Append(":CanPt"); // 15 Ds: Pt
565  listVar.Append(":CanDCAProngToProng"); // 16 Ds: DCA prong-to-prong (K0s,K)
566  listVar.Append(":CanCosThetaStarK0s"); // 17 Ds: cosine theta* K0s
567  listVar.Append(":CanCosThetaStarBach"); // 18 Ds: cosine theta* K+
568  listVar.Append(":CanCosPA"); // 19 Ds: cosine pointing angle
569  listVar.Append(":CanCosPAxy"); // 20 Ds: cosine pointing angle XY
570  listVar.Append(":CanDLengthXY"); // 21 Ds: decay length XY
571  listVar.Append(":CanNormDLengthXY"); // 22 Ds: normalised decay length XY
572  listVar.Append(":CanDLength3D"); // 23 Ds: decay length 3D
573  listVar.Append(":CanNormDLength3D"); // 24 Ds: normalised decay length 3D
574  listVar.Append(":CanSigmaVtx"); // 25 Ds: sigma vertex
575  listVar.Append(":CanNormTopoBach"); // 26 Ds: difference between measured and expected normalised d0(K+)
576  listVar.Append(":CanNormTopoK0s"); // 27 Ds: difference between measured and expected normalised d0(K0s)
577  if (fReadMC) {
578  listVar.Append(":mcTruth"); // 28 MC: MC informations: 10 (signal), 20 (D+ reflection),
579  // -1 (backg w/ true K0S), -2 (backg w/ false K0S)
580  }
581 
582 
583  fOutputNtuple = new TNtuple(GetOutputSlot(4)->GetContainer()->GetName(), "TNtuple for reconstructed candidates on real data", listVar.Data());
584 
585  }
586 
587 
588 
589  PostData(2, fCounter);
590  PostData(3, fOutputSele);
591  if (!fFillNtuple) {
592  PostData(4, fOutputCand);
593  PostData(5, fOutputPID);
594  if (fReadMC)
595  PostData(6, fOutputMC);
596  } else {
597  PostData(4, fOutputNtuple);
598  }
599 
600  return;
601 }
602 
603 
604 
605 
606 
607 //__________________________________________________________________________
609 {
613  AliDebug(1, "AliAnalysisTaskSEDstoK0sK::UserExec()");
614 
615 
616  //---------------------------------------------------------------------------
617  // - Load the event and the CascadesHF branch
618  //---------------------------------------------------------------------------
619 
620  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
621  fHisNEvents->Fill(0); // all events
622 
623 
624  if(fAODProtection>=0){
625  // Protection against different number of events in the AOD and deltaAOD
626  // In case of discrepancy the event is rejected.
627  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
628  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
629  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
630  fHisNEvents->Fill(2);
631  return;
632  }
633  fHisNEvents->Fill(1);
634  }
635 
636 
637  TClonesArray *arrayCascadesHF = 0;
638  if(!aod && AODEvent() && IsStandardAOD()) {
639  // In case there is an AOD handler writing a standard AOD, use the AOD
640  // event in memory rather than the input (ESD) event.
641  aod = dynamic_cast<AliAODEvent*> (AODEvent());
642  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
643  // have to taken from the AOD event hold by the AliAODExtension
644  AliAODHandler *aodHandler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
645  if (aodHandler->GetExtensions()) {
646  AliAODExtension *ext = (AliAODExtension*) aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
647  AliAODEvent *aodFromExt = ext->GetAOD();
648  arrayCascadesHF = (TClonesArray*) aodFromExt->GetList()->FindObject("CascadesHF");
649  }
650  } else if (aod) {
651  arrayCascadesHF = (TClonesArray*) aod->GetList()->FindObject("CascadesHF");
652  }
653 
654  if (!aod || !arrayCascadesHF) {
655  AliWarning("AliAnalysisTaskSEDstoK0sK::UserExec: CascadesHF branch not found!");
656  return;
657  }
658 
659 
660  // Load MC particles
661  TClonesArray *mcArray = 0;
662  AliAODMCHeader *mcHeader = 0;
663 
664  if (fReadMC) {
665  mcArray = (TClonesArray*) aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
666  if (!mcArray) {
667  AliWarning("AliAnalysisTaskSEDstoK0sK::UserExec: MC particles branch not found!");
668  return;
669  }
670  mcHeader = (AliAODMCHeader*) aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
671  if (!mcHeader) {
672  AliWarning("AliAnalysisTaskSEDstoK0sK::UserExec: MC header branch not found!");
673  return;
674  }
675 
676  // FillMCAcceptanceHistos has to take place here
677  }
678 
679 
680 
681  //---------------------------------------------------------------------------
682  // - Event selection
683  //---------------------------------------------------------------------------
684 
685  // Fix for temporary bug in ESDfilter
686  // The AODs with null vertex pointer didn't pass the PhysSel
687  if (!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
688 
689  fHisNEvents->Fill(3);
691 
692 
693  Float_t nTracks = aod->GetNumberOfTracks();
695 
696  fHisCentrality[0]->Fill(evCentr);
697  fHisCentralityVsMult[0]->Fill(nTracks, evCentr);
698 
699 
700  if ( !(fAnalysisCuts->IsEventSelected(aod)) ) {
701 
702  if (fAnalysisCuts->GetWhyRejection() == 5) fHisNEvents->Fill(4); // Rejection: Trigger
703  if (fAnalysisCuts->GetWhyRejection() == 0) fHisNEvents->Fill(5); // Rejection: Vertex reco
704  if (fAnalysisCuts->GetWhyRejection() == 1) fHisNEvents->Fill(6); // Rejection: Pileup
705  if (fAnalysisCuts->GetWhyRejection() == 2) fHisNEvents->Fill(7); // Rejection: Centrality
706  if (fAnalysisCuts->GetWhyRejection() == 6) fHisNEvents->Fill(8); // Rejection: VtxZ fiducial acc.
707  if (fAnalysisCuts->GetWhyRejection() == 7) fHisNEvents->Fill(9); // Rejection: Physics selection
708 
709  if (fAnalysisCuts->GetWhyRejection() == 2) {
710  fHisCentrality[2]->Fill(evCentr);
711  fHisCentralityVsMult[2]->Fill(nTracks, evCentr);
712  }
713 
714  return;
715  }
716 
717 
718  fHisNEvents->Fill(10);
719  fHisCentrality[1]->Fill(evCentr);
720  fHisCentralityVsMult[1]->Fill(nTracks, evCentr);
721 
722 
723  Int_t nCascades = arrayCascadesHF->GetEntriesFast();
724  AliDebug(1, Form("AliAnalysisTaskSEDstoK0sK::UserExec: %d Cascades found for this event!", nCascades));
725 
726 
727 
728 
729  //---------------------------------------------------------------------------
730  // - Analysis: Loop over the candidates
731  //---------------------------------------------------------------------------
732 
733  // vHF object is needed to call the method that refills the missing info of the candidates
734  // if they have been deleted in dAOD reconstruction phase
735  // in order to reduce the size of the file
737 
738 
739 
740  for (Int_t iCasc=0; iCasc<nCascades; iCasc++)
741  {
742 
743  //-----------------------------------------
744  // - Get candidates and fill missing info
745  //-----------------------------------------
746  AliAODRecoCascadeHF *dCan = (AliAODRecoCascadeHF*) arrayCascadesHF->UncheckedAt(iCasc);
747 
748  if (fUseSelectionBit) {
749  if ( !(dCan->CheckCascadeFlags(AliRDHFCuts::kDstoK0sCuts)) ) continue;
750  }
751 
752  fHisNEvents->Fill(11);
753 
754 
755  // Fill the data members of the candidate only if they are empty.
756  if (!(vHF->FillRecoCasc(aod, dCan, kFALSE, kTRUE))) {
757  fHisNEvents->Fill(12);
758  continue;
759  }
760 
761 
762  // Reconstruct secondary vertex if it does not already exist (standard dAOD)
763  if (dCan->GetIsFilled()==1) {
764  if (!(vHF->RecoSecondaryVertexForCascades(aod, dCan))) {
765  fHisNEvents->Fill(13);
766  continue;
767  }
768  }
769 
770  fHisNEvents->Fill(14);
771 
772 
773 
774  //-----------------------------------------
775  // - Select candidates step by step
776  //-----------------------------------------
777  Double_t ptCand = dCan->Pt();
778  Double_t rapidity = dCan->Y(431);
779 
780  fHisRapidity->Fill(ptCand, rapidity);
781 
782 
783  if (!(fAnalysisCuts->IsInFiducialAcceptance(ptCand, rapidity))) continue;
784  fHisNEvents->Fill(15);
785 
786 
787  if (!(fAnalysisCuts->IsSelected(dCan, AliRDHFCuts::kTracks, aod))) continue;
788  fHisNEvents->Fill(16);
789 
790  // FillTheTree(dCan, mcArray);
791 
792 
793  if (!(fAnalysisCuts->IsSelected(dCan, AliRDHFCuts::kCandidate, aod))) continue;
794  fHisNEvents->Fill(17);
795 
798 
799 
800  if (!(fAnalysisCuts->IsSelected(dCan, AliRDHFCuts::kPID, aod))) continue;
801  fHisNEvents->Fill(18);
802 
803  fHisRapiditySel->Fill(ptCand, rapidity);
804 
805  FillHistogramsVar(dCan, AliRDHFCuts::kPID, mcArray);
806  FillHistogramsPID(dCan, AliRDHFCuts::kPID, mcArray);
807 
808 
809  FillTheTree(dCan, aod, mcArray);
810 
811 
812  } // end loop cascades
813 
814 
815 
816  delete vHF;
817 
818 
819  PostData(2, fCounter);
820  PostData(3, fOutputSele);
821  if (!fFillNtuple) {
822  PostData(4, fOutputCand);
823  PostData(5, fOutputPID);
824  if (fReadMC)
825  PostData(6, fOutputMC);
826  } else {
827  PostData(4, fOutputNtuple);
828  }
829 
830  return;
831 }
832 
833 
834 
835 
836 //__________________________________________________________________________
838 {
842  AliDebug(1, "AliAnalysisTaskSEDstoK0sK::Terminate()");
843 
844 
845  // Output slot #3
846  fOutputSele = dynamic_cast<TList*> (GetOutputData(3));
847  if (!fOutputSele) {
848  AliWarning("ERROR: fOutputSele not available");
849  return;
850  }
851 
852 
853  if (!fFillNtuple) {
854 
855  // Output slot #4
856  fOutputCand = dynamic_cast<TList*> (GetOutputData(4));
857  if (!fOutputCand) {
858  AliWarning("ERROR: fOutputCand not available");
859  return;
860  }
861 
862  // Output slot #5
863  fOutputPID = dynamic_cast<TList*> (GetOutputData(5));
864  if (!fOutputPID) {
865  AliWarning("ERROR: fOutputPID not available");
866  return;
867  }
868 
869  if (fReadMC) {
870  // Output slot #6
871  fOutputMC = dynamic_cast<TList*> (GetOutputData(6));
872  if (!fOutputMC) {
873  AliWarning("ERROR: fOutputMC not available");
874  return;
875  }
876  }
877 
878  } else {
879 
880  // Output slot #4
881  fOutputNtuple = dynamic_cast<TNtuple*> (GetOutputData(4));
882  if (!fOutputNtuple) {
883  AliWarning("ERROR: fOutputNtuple not available");
884  return;
885  }
886 
887  }
888 
889  return;
890 }
891 
892 
893 
894 
895 //__________________________________________________________________________
897 {
901 
902 
903  if (nBins<=kMaxPtBins) {
904 
905  fNPtBins = nBins;
906  for (Int_t ipt=0; ipt<fNPtBins+1; ipt++) fPtLimits[ipt] = limitsPt[ipt];
907  for (Int_t ipt=fNPtBins+1; ipt<kMaxPtBins+1; ipt++) fPtLimits[ipt] = 99999999.;
908 
909  } else {
910 
911  // More Pt bins than allowed by the class header
912  AliWarning(Form("Maximum number of Pt bins reached: %d", kMaxPtBins));
914  fPtLimits[0] = 2.;
915  fPtLimits[1] = 4.;
916  fPtLimits[2] = 6.;
917  fPtLimits[3] = 8.;
918  for (Int_t ipt=4; ipt<kMaxPtBins+1; ipt++) fPtLimits[ipt] = 99999999.;
919 
920  }
921 
922  if (fDebug>1) {
923  AliInfo(Form("Number of Pt bins for the analysis: %d", fNPtBins));
924  for (Int_t ipt=0; ipt<fNPtBins; ipt++) {
925  AliInfo(Form(" Pt bin %02d: [%8.1f, %8.1f]", ipt, fPtLimits[ipt], fPtLimits[ipt+1]));
926  }
927  }
928 
929  return;
930 }
931 
932 
933 
934 
935 //__________________________________________________________________________
937 {
941 
942 
943  if (nCuts==kNTupleVars) {
944 
945  for (Int_t ivar=0; ivar<kNTupleVars; ivar++) {
946  fCutsMinTupleVars[ivar] = minCuts[ivar];
947  fCutsMaxTupleVars[ivar] = maxCuts[ivar];
948  }
949 
950  } else {
951 
952  // Different number of cuts than the number of variables in the tuple
953  AliWarning(Form("%d cuts are set for %d variables in the tuple: default initialisation", nCuts, kNTupleVars));
954 
955  for (Int_t ivar=0; ivar<kNTupleVars; ivar++) {
956  fCutsMinTupleVars[ivar] = -999.;
957  fCutsMaxTupleVars[ivar] = 999.;
958  }
959 
960  }
961 
962 
963  return;
964 }
965 
966 
967 
968 
969 //__________________________________________________________________________
971 {
982  AliDebug(1, "AliAnalysisTaskSEDstoK0sK::FillHistogramsVar()");
983 
984 
985  if (fFillNtuple) return;
986 
987 
988  //---------------------------------------------------------------------------
989  // - Monte Carlo: check the candidate at MC level
990  //---------------------------------------------------------------------------
991  Bool_t okMCsignal = kFALSE;
992  Bool_t okMCbackg = kFALSE;
993  Bool_t okMCreflec = kFALSE;
994 
995  if (selFlag==AliRDHFCuts::kPID && fReadMC) {
996  if (!mcArray) return;
997 
998  if (MatchToMCDstoK0sKSignal(dCan, mcArray)>=0) {
999  AliDebug(2, "This candidate matches a MC signal of Ds+ -> K0S + K+");
1000  okMCsignal = kTRUE;
1001  } else if (MatchToMCDplustoK0spiSignal(dCan, mcArray)>=0) {
1002  AliDebug(2, "This candidate matches a MC signal of D+ -> K0S + pi+");
1003  okMCreflec = kTRUE;
1004  } else {
1005  AliDebug(2, "This candidate is a MC background");
1006  okMCbackg = kTRUE;
1007  }
1008  }
1009 
1010 
1011 
1012  //---------------------------------------------------------------------------
1013  // - Fill the histograms depending on the selection level
1014  // - if 'kCandidate' level, fill only index=0 (histograms 'kCand')
1015  // - if 'kPID' level and not reading MC, fill only index=1 (histograms 'kPID')
1016  // - if 'kPID' level and reading MC, fill index=1 ('kPID'), 2 ('mcSig'), 3 ('mcBackg'), 4 ('mcDplusRefl')
1017  //---------------------------------------------------------------------------
1018  Float_t ptCand = dCan->Pt();
1019 
1020 
1021  for (Int_t index=0; index<5; index++) {
1022 
1023  if (index==0) {
1024  if (selFlag != AliRDHFCuts::kCandidate) continue;
1025  } else if (index==1) {
1026  if (selFlag != AliRDHFCuts::kPID) continue;
1027  } else if (index==2) {
1028  if (!okMCsignal) continue;
1029  } else if (index==3) {
1030  if (!okMCbackg) continue;
1031  } else if (index==4) {
1032  if (!okMCreflec) continue;
1033  }
1034 
1035  fHisInvMassDs[index] ->Fill( ptCand, dCan->InvMassDstoK0sK() );
1036  fHisInvMassDplus[index] ->Fill( ptCand, dCan->InvMassDplustoK0spi() );
1037  fHisInvMassK0s[index] ->Fill( ptCand, dynamic_cast<AliAODv0*>(dCan->Getv0())->MassK0Short());
1038  fHisPtK0s[index] ->Fill( ptCand, dCan->PtProng(1) );
1039  fHisPtBachelor[index] ->Fill( ptCand, dCan->PtProng(0) );
1040  fHisImpParK0s[index] ->Fill( ptCand, dCan->Getd0Prong(1) );
1041  fHisImpParBach[index] ->Fill( ptCand, dCan->Getd0Prong(0) );
1042  fHisCTauK0s[index] ->Fill( ptCand, dynamic_cast<AliAODv0*>(dCan->Getv0())->Ct(310, dCan->GetOwnPrimaryVtx()));
1043  fHisCosPointingDs[index] ->Fill( ptCand, dCan->CosPointingAngle() );
1044  fHisCosPointingXYDs[index] ->Fill( ptCand, dCan->CosPointingAngleXY() );
1045  fHisCosThetaStarK0s[index] ->Fill( ptCand, dCan->CosThetaStar(1, 431, 321, 310) );
1046  fHisCosThetaStarBach[index] ->Fill( ptCand, dCan->CosThetaStar(0, 431, 321, 310) );
1047  fHisDCAK0sBach[index] ->Fill( ptCand, dCan->GetDCA() );
1048  fHisDecayLxyDs[index] ->Fill( ptCand, dCan->DecayLengthXY() );
1049  fHisNormDecayLxyDs[index] ->Fill( ptCand, dCan->NormalizedDecayLengthXY() );
1050 
1051  }
1052 
1053 
1054  return;
1055 }
1056 
1057 
1058 
1059 
1060 //__________________________________________________________________________
1062 {
1066  AliDebug(2, "AliAnalysisTaskSEDstoK0sK::FillHistogramsPID()");
1067 
1068 
1069  Double_t nSigmaTPC = -999;
1070  Double_t nSigmaTOF = -999.;
1071  Float_t ptCand = dCan->Pt();
1072 
1073  if (fAnalysisCuts->GetIsUsePID()) {
1075  if (pidhf->CheckTPCPIDStatus(dCan->GetBachelor())) pidhf->GetnSigmaTPC(dCan->GetBachelor(), AliPID::kKaon, nSigmaTPC);
1076  if (pidhf->CheckTOFPIDStatus(dCan->GetBachelor())) pidhf->GetnSigmaTOF(dCan->GetBachelor(), AliPID::kKaon, nSigmaTOF);
1077  }
1078 
1079 
1080  if (selFlag==AliRDHFCuts::kCandidate) {
1081 
1082  //---------------------------------------------------------------------------
1083  // - Topological selection levels
1084  //---------------------------------------------------------------------------
1085  fHisPidTPCKaonVsPt->Fill(ptCand, nSigmaTPC);
1086  fHisPidTOFKaonVsPt->Fill(ptCand, nSigmaTOF);
1087  fHisPidTPCTOFKaon->Fill(nSigmaTPC, nSigmaTOF);
1088 
1089  } else if (selFlag==AliRDHFCuts::kPID || selFlag==AliRDHFCuts::kAll) {
1090 
1091  //---------------------------------------------------------------------------
1092  // - Topological+PID selection levels
1093  //---------------------------------------------------------------------------
1094  fHisPidTPCKaonVsPtSel->Fill(ptCand, nSigmaTPC);
1095  fHisPidTOFKaonVsPtSel->Fill(ptCand, nSigmaTOF);
1096  fHisPidTPCTOFKaonSel->Fill(nSigmaTPC, nSigmaTOF);
1097 
1098  //-----------------------------------------
1099  // - Check the candidate at MC level
1100  //-----------------------------------------
1101  if (fReadMC && mcArray) {
1102  if (MatchToMCDstoK0sKSignal(dCan, mcArray)>=0) fHisNEvents->Fill(19);
1103  else if (MatchToMCDplustoK0spiSignal(dCan, mcArray)>=0) fHisNEvents->Fill(20);
1104  }
1105 
1106  }
1107 
1108 
1109  return;
1110 }
1111 
1112 
1113 
1114 
1115 //__________________________________________________________________________
1117 {
1121  AliDebug(2, "AliAnalysisTaskSEDstoK0sK::FillTheTree()");
1122 
1123 
1124  if (!fFillNtuple) return;
1125 
1126  AliAODv0 *v0 = (AliAODv0*) dCan->Getv0();
1127  if (!v0) return;
1128 
1129 
1130  //---------------------------------------------------------------------------
1131  // - Get PID informations if possible
1132  //---------------------------------------------------------------------------
1133  Double_t nSigmaTPCpion = -999;
1134  Double_t nSigmaTPCkaon = -999;
1135  Double_t nSigmaTPCproton = -999;
1136  Double_t nSigmaTOFpion = -999;
1137  Double_t nSigmaTOFkaon = -999;
1138  Double_t nSigmaTOFproton = -999;
1139  if (fAnalysisCuts->GetIsUsePID()) {
1141  fAnalysisCuts->GetPidHF()->GetnSigmaTPC(dCan->GetBachelor(), AliPID::kPion, nSigmaTPCpion);
1142  fAnalysisCuts->GetPidHF()->GetnSigmaTPC(dCan->GetBachelor(), AliPID::kKaon, nSigmaTPCkaon);
1143  fAnalysisCuts->GetPidHF()->GetnSigmaTPC(dCan->GetBachelor(), AliPID::kProton, nSigmaTPCproton);
1144  }
1146  fAnalysisCuts->GetPidHF()->GetnSigmaTOF(dCan->GetBachelor(), AliPID::kPion, nSigmaTOFpion);
1147  fAnalysisCuts->GetPidHF()->GetnSigmaTOF(dCan->GetBachelor(), AliPID::kKaon, nSigmaTOFkaon);
1148  fAnalysisCuts->GetPidHF()->GetnSigmaTOF(dCan->GetBachelor(), AliPID::kProton, nSigmaTOFproton);
1149  }
1150  }
1151 
1152 
1153  //---------------------------------------------------------------------------
1154  // - Apply cuts to the tuple variables
1155  //---------------------------------------------------------------------------
1156  Double_t diffIP[2] = {0.}, errdiffIP[2] = {1.};
1157  dCan->Getd0MeasMinusExpProng(0, aod->GetMagneticField(), diffIP[0], errdiffIP[0]);
1158  dCan->Getd0MeasMinusExpProng(1, aod->GetMagneticField(), diffIP[1], errdiffIP[1]);
1159 
1160 
1161  const Int_t nVarTuple = fReadMC ? kNTupleVarsMC : kNTupleVars;
1162 
1163  Float_t variableNtuple[nVarTuple];
1164  variableNtuple[ 0] = v0->MassK0Short(); // "K0sInvMass"
1165  variableNtuple[ 1] = dCan->PtProng(1); // "K0sPt"
1166  variableNtuple[ 2] = v0->RadiusV0(); // "K0sRxy"
1167  variableNtuple[ 3] = v0->Ct(310, dCan->GetOwnPrimaryVtx()); // "K0sCTau"
1168  variableNtuple[ 4] = dCan->CosV0PointingAngle(); // "K0sCosPA"
1169  variableNtuple[ 5] = dCan->Getd0Prong(1); // "K0sd0"
1170  variableNtuple[ 6] = v0->DcaPosToPrimVertex(); // "K0sd0DauPos"
1171  variableNtuple[ 7] = v0->DcaNegToPrimVertex(); // "K0sd0DauNeg"
1172  variableNtuple[ 8] = v0->DcaV0Daughters(); // "K0sDCAdau"
1173  variableNtuple[ 9] = v0->PtProng(0); // "K0sPtDauPos"
1174  variableNtuple[10] = v0->PtProng(1); // "K0sPtDauNeg"
1175  variableNtuple[11] = dCan->PtProng(0); // "BachPt"
1176  variableNtuple[12] = dCan->Getd0Prong(0); // "Bachd0"
1177  variableNtuple[13] = dCan->InvMassDstoK0sK(); // "CanInvMassDs"
1178  variableNtuple[14] = dCan->InvMassDplustoK0spi(); // "CanInvMassDplus"
1179  variableNtuple[15] = dCan->Pt(); // "CanPt"
1180  variableNtuple[16] = dCan->GetDCA(); // "CanDCAProngToProng"
1181  variableNtuple[17] = dCan->CosThetaStar(1, 431, 321, 310); // "CanCosThetaStarK0s"
1182  variableNtuple[18] = dCan->CosThetaStar(0, 431, 321, 310); // "CanCosThetaStarBach"
1183  variableNtuple[19] = dCan->CosPointingAngle(); // "CanCosPA"
1184  variableNtuple[20] = dCan->CosPointingAngleXY(); // "CanCosPAxy"
1185  variableNtuple[21] = dCan->DecayLengthXY(); // "CanDLengthXY"
1186  variableNtuple[22] = dCan->NormalizedDecayLengthXY(); // "CanNormDLengthXY"
1187  variableNtuple[23] = dCan->DecayLength(); // "CanDLength3D"
1188  variableNtuple[24] = dCan->NormalizedDecayLength(); // "CanNormDLength3D"
1189  variableNtuple[25] = ComputeSigmaVert(aod, dCan); // "CanSigmaVtx"
1190  variableNtuple[26] = diffIP[0]/errdiffIP[0]; // "CanNormTopoBach"
1191  variableNtuple[27] = diffIP[1]/errdiffIP[1]; // "CanNormTopoK0s"
1192 
1193 
1194  for (Int_t ivar=0; ivar<kNTupleVars; ivar++) {
1195  if ( (variableNtuple[ivar]<fCutsMinTupleVars[ivar]) || (variableNtuple[ivar]>fCutsMaxTupleVars[ivar]) ) {
1196  AliDebug(2, "One of the tuple variable does not pass the selection cut");
1197  return;
1198  }
1199  }
1200 
1201 
1202  //---------------------------------------------------------------------------
1203  // - Monte Carlo: check the candidate at MC level
1204  //---------------------------------------------------------------------------
1205  if (fReadMC) {
1206  Bool_t okMCsignal = (MatchToMCDstoK0sKSignal(dCan, mcArray)>=0) ? kTRUE : kFALSE;
1207  Bool_t okMCreflec = (MatchToMCDplustoK0spiSignal(dCan, mcArray)>=0) ? kTRUE : kFALSE;
1208 
1209  variableNtuple[28] = okMCsignal ? 10
1210  : okMCreflec ? 20
1211  : -2;
1212 
1213  if (!okMCsignal && !okMCreflec) {
1214  AliDebug(2, "The candidate is a combinatorial background, check if the V0 is a true K0S");
1215  Int_t pdgDgK0stoPions[2] = {211, 211};
1216  if (v0->MatchToMC(310, mcArray, 2, pdgDgK0stoPions)>=0) variableNtuple[28]++;
1217  }
1218 
1219  } // end if fReadMC
1220 
1221 
1222  //---------------------------------------------------------------------------
1223  // - Fill the tuple
1224  //---------------------------------------------------------------------------
1225  fOutputNtuple->Fill(variableNtuple);
1226 
1227 
1228  return;
1229 }
1230 
1231 
1232 
1233 
1234 //__________________________________________________________________________
1236 {
1241 
1242 
1243  // - Get bachelor and V0 from the cascade
1244  AliESDtrack *esdB = new AliESDtrack(dCan->GetBachelor());
1245 
1246  AliNeutralTrackParam *trackV0 = 0;
1247 
1248  const AliVTrack *trackVV0 = dynamic_cast<const AliVTrack*>(dCan->Getv0());
1249  if (trackVV0) {
1250  trackV0 = new AliNeutralTrackParam(trackVV0);
1251  }
1252 
1253 
1254  if (!esdB || !trackV0) {
1255  if (esdB) delete esdB;
1256  if (trackV0) delete trackV0;
1257  return -999.;
1258  }
1259 
1260 
1261  // - Build ESD primary vertex
1262  AliVertexerTracks vertexer(aod->GetMagneticField());
1263 
1264  Double_t pos[3], cov[6];
1265  AliAODVertex *aodV = (AliAODVertex*) aod->GetPrimaryVertex();
1266  aodV->GetXYZ(pos);
1267  aodV->GetCovarianceMatrix(cov);
1268  Double_t chi2 = aodV->GetChi2();
1269  Int_t nContr = aodV->GetNContributors();
1270 
1271  AliESDVertex vprim(pos, cov, chi2, nContr);
1272  vertexer.SetVtxStart(&vprim);
1273 
1274 
1275  // - Assign TObjArray of daughters
1276  TObjArray twoTrackArray(2);
1277  twoTrackArray.AddAt(esdB, 0);
1278  twoTrackArray.AddAt(trackV0, 1);
1279 
1280 
1281  // - Build ESD secondary vertex and get dispersion
1282  AliESDVertex *secVert = vertexer.VertexForSelectedESDTracks(&twoTrackArray, kFALSE, kTRUE, kFALSE);
1283  Double_t disp = secVert->GetDispersion();
1284 
1285 
1286  twoTrackArray.Delete();
1287  delete secVert;
1288 
1289  return disp;
1290 }
1291 
1292 
1293 
1294 
1295 //__________________________________________________________________________
1297 {
1302 
1303 
1304  // - Get bachelor and V0 from the cascade
1305  AliAODv0 *v0 = dynamic_cast<AliAODv0*> (dCan->Getv0());
1306  AliAODTrack *bach = dynamic_cast<AliAODTrack*> (dCan->GetBachelor());
1307  if (!v0 || !bach) {
1308  return -999;
1309  }
1310 
1311 
1312  // - Get cascade informations necessary to compute the rest frame
1313  Double_t eCasc = dCan->EProng(0, 321) + dCan->EProng(1, 310);
1314  Double_t pxCasc = dCan->PxProng(0) + dCan->PxProng(1);
1315  Double_t pyCasc = dCan->PyProng(0) + dCan->PyProng(1);
1316  Double_t pzCasc = dCan->PzProng(0) + dCan->PzProng(1);
1317  Double_t bxCasc = pxCasc/eCasc;
1318  Double_t byCasc = pyCasc/eCasc;
1319  Double_t bzCasc = pzCasc/eCasc;
1320 
1321 
1322  // - Boost the K+ and K0s in the cascade rest frame
1323  TVector3 vecK0sCascFrame;
1324  TLorentzVector vecK0s(dCan->PxProng(1), dCan->PyProng(1), dCan->PzProng(1), dCan->EProng(1, 310));
1325  vecK0s.Boost(-bxCasc, -byCasc, -bzCasc);
1326  vecK0s.Boost(vecK0sCascFrame);
1327  vecK0sCascFrame = vecK0s.BoostVector();
1328 
1329  TVector3 vecBachCascFrame;
1330  TLorentzVector vecBach(dCan->PxProng(0), dCan->PyProng(0), dCan->PzProng(0), dCan->EProng(0, 321));
1331  vecBach.Boost(-bxCasc, -byCasc, -bzCasc);
1332  vecBach.Boost(vecBachCascFrame);
1333  vecBachCascFrame = vecBach.BoostVector();
1334 
1335 
1336  // - Compute the cosine of the angle between K+ and K0s
1337  Double_t scalProdBachK0s = vecBachCascFrame.Dot(vecK0sCascFrame);
1338  Double_t normBach = vecBachCascFrame.Mag();
1339  Double_t normK0s = vecK0sCascFrame.Mag();
1340  Double_t cosAngleCascFrame = scalProdBachK0s/(normBach*normK0s);
1341 
1342 
1343  return cosAngleCascFrame;
1344 }
1345 
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:315
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:247
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:258
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
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:248
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:270
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:249
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:280
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