AliPhysics  32b88a8 (32b88a8)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskSEDmesonsFilterCJ.cxx
Go to the documentation of this file.
1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 //
16 //
17 // Task for selecting D mesons to be used as an input for D-Jet correlations
18 //
19 //-----------------------------------------------------------------------
20 // Authors:
21 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
22 // A. Grelli (Utrecht University) a.grelli@uu.nl
23 // X. Zhang (LBNL) XMZhang@lbl.gov
24 // S. Aiola (Yale University) salvatore.aiola@cern.ch
25 //-----------------------------------------------------------------------
26 
27 #include <TDatabasePDG.h>
28 #include <TParticle.h>
29 #include <TVector3.h>
30 #include <TROOT.h>
31 #include <TH3F.h>
32 
33 #include "AliLog.h"
35 #include "AliRDHFCutsD0toKpi.h"
36 #include "AliAODMCHeader.h"
37 #include "AliAODHandler.h"
38 #include "AliAnalysisManager.h"
39 #include "AliLog.h"
40 #include "AliAODVertex.h"
41 #include "AliAODRecoCascadeHF.h"
43 #include "AliESDtrack.h"
44 #include "AliAODMCParticle.h"
45 #include "AliEmcalParticle.h"
46 #include "AliParticleContainer.h"
47 #include "AliAnalysisDataSlot.h"
48 #include "AliAnalysisDataContainer.h"
49 #include "AliStack.h"
50 
52 
54 
55 //_______________________________________________________________________________
58  fUseMCInfo(kFALSE),
59  fUseReco(kTRUE),
60  fCandidateType(0),
61  fCandidateName(""),
62  fPDGmother(0),
63  fNProngs(0),
64  fBranchName(""),
65  fCuts(0),
66  fMinMass(0.),
67  fMaxMass(0.),
68  fInhibitTask(kFALSE),
69  fCombineDmesons(kFALSE),
70  fMultCand(kFALSE),
71  fAnalyseCand(0),
72  fRejectQuarkNotFound(kTRUE),
73  fRejectDfromB(kTRUE),
74  fKeepOnlyDfromB(kFALSE),
75  fAodEvent(0),
76  fArrayDStartoD0pi(0),
77  fMCarray(0),
78  fCandidateArray(0),
79  fSideBandArray(0),
80  fCombinedDmesons(0),
81  fCombinedDmesonsBkg(0),
82  fNCand(0),
83  fNSBCand(0),
84  fHistStat(0),
85  fHistNSBCandEv(0),
86  fHistNCandEv(0),
87  fHistImpParS(0),
88  fHistImpParB(0),
89  fHistPtPion(0),
90  fHistInvMassPtD(0),
91  fHistInvMassS(0),
92  fHistInvMassB(0),
93  fHistAlphaDDS(0),
94  fHistAlphaDpisS(0),
95  fHistAlphaDpiS(0),
96  fHistAlphaDKS(0),
97  fHistAlphaDDB(0),
98  fHistAlphaDpisB(0),
99  fHistAlphaDpiB(0),
100  fHistAlphaDKB(0),
101  fHistDeltaRDDS(0),
102  fHistDeltaRDpisS(0),
103  fHistDeltaRDpiS(0),
104  fHistDeltaRDKS(0),
105  fHistDeltaRDDB(0),
106  fHistDeltaRDpisB(0),
107  fHistDeltaRDpiB(0),
108  fHistDeltaRDKB(0),
109  fHistAlphaDpiR(0),
110  fHistAlphaDKR(0),
111  fHistDeltaRDpiR(0),
112  fHistDeltaRDKR(0)
113 {
114  //
115  // Default constructor
116  //
117 
118  for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
119  for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
120 
121  fNeedEmcalGeom = kFALSE;
122 }
123 
124 //_______________________________________________________________________________
126  AliAnalysisTaskEmcal(name, kTRUE),
127  fUseMCInfo(kFALSE),
128  fUseReco(kTRUE),
129  fCandidateType(candtype),
130  fCandidateName(""),
131  fPDGmother(0),
132  fNProngs(0),
133  fBranchName(""),
134  fCuts(cuts),
135  fMinMass(0.),
136  fMaxMass(0.),
137  fInhibitTask(kFALSE),
138  fCombineDmesons(kFALSE),
139  fMultCand(kFALSE),
140  fAnalyseCand(0),
141  fRejectQuarkNotFound(kTRUE),
142  fRejectDfromB(kTRUE),
143  fKeepOnlyDfromB(kFALSE),
144  fAodEvent(0),
145  fArrayDStartoD0pi(0),
146  fMCarray(0),
147  fCandidateArray(0),
148  fSideBandArray(0),
149  fCombinedDmesons(0),
150  fCombinedDmesonsBkg(0),
151  fNCand(0),
152  fNSBCand(0),
153  fHistStat(0),
154  fHistNSBCandEv(0),
155  fHistNCandEv(0),
156  fHistImpParS(0),
157  fHistImpParB(0),
158  fHistPtPion(0),
159  fHistInvMassPtD(0),
160  fHistInvMassS(0),
161  fHistInvMassB(0),
162  fHistAlphaDDS(0),
163  fHistAlphaDpisS(0),
164  fHistAlphaDpiS(0),
165  fHistAlphaDKS(0),
166  fHistAlphaDDB(0),
167  fHistAlphaDpisB(0),
168  fHistAlphaDpiB(0),
169  fHistAlphaDKB(0),
170  fHistDeltaRDDS(0),
171  fHistDeltaRDpisS(0),
172  fHistDeltaRDpiS(0),
173  fHistDeltaRDKS(0),
174  fHistDeltaRDDB(0),
175  fHistDeltaRDpisB(0),
176  fHistDeltaRDpiB(0),
177  fHistDeltaRDKB(0),
178  fHistAlphaDpiR(0),
179  fHistAlphaDKR(0),
180  fHistDeltaRDpiR(0),
181  fHistDeltaRDKR(0)
182 {
183  //
184  // Constructor. Initialization of Inputs and Outputs
185  //
186 
187  Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
188 
189  fNeedEmcalGeom = kFALSE;
190 
191  for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
192  for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
193 
194  const Int_t nptbins = fCuts->GetNPtBins();
195  Float_t defaultSigmaD013[13] = { 0.012, 0.012, 0.012, 0.015, 0.015, 0.018, 0.018, 0.020, 0.020, 0.030, 0.030, 0.037, 0.040 };
196 
197  switch (fCandidateType) {
198  case 0 :
199  fCandidateName = "D0";
200  fPDGmother = 421;
201  fNProngs = 2;
202  fPDGdaughters[0] = 211; // pi
203  fPDGdaughters[1] = 321; // K
204  fPDGdaughters[2] = 0; // empty
205  fPDGdaughters[3] = 0; // empty
206  fBranchName = "D0toKpi";
207  break;
208  case 1 :
209  fCandidateName = "DStar";
210  fPDGmother = 413;
211  fNProngs = 3;
212  fPDGdaughters[1] = 211; // pi soft
213  fPDGdaughters[0] = 421; // D0
214  fPDGdaughters[2] = 211; // pi fromD0
215  fPDGdaughters[3] = 321; // K from D0
216  fBranchName = "Dstar";
217 
218  if (nptbins<=13) {
219  for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
220  }
221  else {
222  AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
223  }
224  break;
225  default :
226  Warning("AliAnalysisTaskSEDmesonsFilterCJ", "%d not accepted!!", fCandidateType);
227  break;
228  }
229 
232 
233  DefineOutput(1, TList::Class()); // histos
234  DefineOutput(2, AliRDHFCuts::Class()); // my cuts
235  DefineOutput(3, TClonesArray::Class()); //array of candidates
236  DefineOutput(4, TClonesArray::Class()); //array of SB candidates
237  DefineOutput(5, TClonesArray::Class()); //array of candidates and event tracks
238  DefineOutput(6, TClonesArray::Class()); //array of SB candidates and event tracks
239 }
240 
241 //_______________________________________________________________________________
243 {
244  //
245  // Destructor
246  //
247 
248  Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");
249 
250  if (fCuts) {
251  delete fCuts;
252  fCuts = 0;
253  }
254 
255  if (fCandidateArray) {
256  delete fCandidateArray;
257  fCandidateArray = 0;
258  }
259 
260  if (fSideBandArray) {
261  delete fSideBandArray;
262  fSideBandArray = 0;
263  }
264 
265  if (fCombinedDmesons) {
266  delete fCombinedDmesons;
267  fCombinedDmesons = 0;
268  }
269 
270  if (fCombinedDmesonsBkg) {
271  delete fCombinedDmesonsBkg;
273  }
274 }
275 
276 //_______________________________________________________________________________
278 {
279  //
280  // Initialization
281  //
282 
283  Info("AnalysisTaskSEDmesonsForJetCorrelations::Init()", "Entering method");
284 
285  switch (fCandidateType) {
286  case 0:
287  {
288  AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
289  copyfCutsDzero->SetName("AnalysisCutsDzero");
290  PostData(2, copyfCutsDzero); // Post the data
291  } break;
292  case 1:
293  {
294  AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
295  copyfCutsDstar->SetName("AnalysisCutsDStar");
296  PostData(2, copyfCutsDstar); // Post the cuts
297  } break;
298  default:
299  return;
300  }
301 
302  return;
303 }
304 
305 //_______________________________________________________________________________
307 {
308  //
309  // Create output objects
310  //
311 
312  Info("UserCreateOutputObjects","CreateOutputObjects of task %s", GetName());
313 
315 
316  DefineHistoForAnalysis(); // define histograms
317 
318  if (fUseReco) {
319  if (fCandidateType == kD0toKpi){
320  fCandidateArray = new TClonesArray("AliAODRecoDecayHF2Prong",10);
321  fSideBandArray = new TClonesArray("AliAODRecoDecayHF2Prong",10);
322  }
323  else if (fCandidateType == kDstartoKpipi) {
324  fCandidateArray = new TClonesArray("AliAODRecoCascadeHF",10);
325  fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",10);
326  }
327  else {
328  AliWarning(Form("Candidate type %d not recognized!", fCandidateType));
329  return;
330  }
331  }
332  else {
333  fCandidateArray = new TClonesArray("AliAODMCParticle",10);
334  fSideBandArray = new TClonesArray("TObject",0); // not used
335  }
336 
337  if (fCombineDmesons) {
338  fCombinedDmesons = new TClonesArray("AliEmcalParticle",50);
339  fCombinedDmesonsBkg = new TClonesArray("AliEmcalParticle",50);
340  }
341  else {
342  fCombinedDmesons = new TClonesArray("TObject",0); // not used
343  fCombinedDmesonsBkg = new TClonesArray("TObject",0); // not used
344  }
345 
346  fCandidateArray->SetOwner();
347  fCandidateArray->SetName(GetOutputSlot(3)->GetContainer()->GetName());
348 
349  //this is used for the DStar side bands and MC!
350  fSideBandArray->SetOwner();
351  fSideBandArray->SetName(GetOutputSlot(4)->GetContainer()->GetName());
352 
353  fCombinedDmesons->SetOwner();
354  fCombinedDmesons->SetName(GetOutputSlot(5)->GetContainer()->GetName());
355 
356  fCombinedDmesonsBkg->SetOwner();
357  fCombinedDmesonsBkg->SetName(GetOutputSlot(6)->GetContainer()->GetName());
358 
359  PostData(1, fOutput);
360  PostData(3, fCandidateArray);
361  PostData(4, fSideBandArray);
362  PostData(5, fCombinedDmesons);
363  PostData(6, fCombinedDmesonsBkg);
364 
365  Info("UserCreateOutputObjects","Data posted for task %s", GetName());
366 }
367 
368 //_______________________________________________________________________________
370 {
371  //
372  // To be executed only once, for the first event
373  //
374 
375  AliDebug(2, "Entering ExecOnce()");
376 
377  if (fInhibitTask) return;
378 
379  // Load the event
380  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
381 
382  if (fAodEvent) {
383  fArrayDStartoD0pi = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(fBranchName.Data()));
384  }
385  else {
386  if (AODEvent() && IsStandardAOD()) {
387 
388  // In case there is an AOD handler writing a standard AOD, use the AOD
389  // event in memory rather than the input (ESD) event.
390  fAodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
391 
392  // in this case the branches in the deltaAOD (AliAOD.VertexingHF.root)
393  // have to taken from the AOD event hold by the AliAODExtension
394  AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
395  if(aodHandler->GetExtensions()) {
396  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
397  AliAODEvent *aodFromExt = ext->GetAOD();
398  fArrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
399  }
400  else {
401  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
402  fInhibitTask = kTRUE;
403  return;
404  }
405  }
406  }
407 
408  if (fArrayDStartoD0pi) {
409  TString objname(fArrayDStartoD0pi->GetClass()->GetName());
410  TClass cls(objname);
411  if (!cls.InheritsFrom("AliAODRecoDecayHF2Prong")) {
412  AliError(Form("%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
413  GetName(), cls.GetName(), fArrayDStartoD0pi->GetName()));
414  fInhibitTask = kTRUE;
415  fArrayDStartoD0pi = 0;
416  return;
417  }
418  }
419  else {
420  AliError(Form("Could not find array %s, skipping the event. Task '%s' will be disabled!", fBranchName.Data(), GetName()));
421  fInhibitTask = kTRUE;
422  return;
423  }
424 
425  if (fUseMCInfo) {
426  fMCarray = dynamic_cast<TClonesArray*>(fAodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
427  if (!fMCarray) {
428  AliError(Form("MC particles not found! Task '%s' will be disabled!", GetName()));
429  fInhibitTask = kTRUE;
430  return;
431  }
432  }
433 
434  if (fCombineDmesons) {
437  }
438 
440 }
441 
442 //_______________________________________________________________________________
444 {
445  //
446  // Analysis execution
447  //
448 
449  if (fInhibitTask) return kFALSE;
450 
451  AliDebug(2, "Entering Run()");
452 
453  //clear the TClonesArray from the previous event
454  fCandidateArray->Clear();
455  fSideBandArray->Clear();
456  fCombinedDmesons->Clear();
457  fCombinedDmesonsBkg->Clear();
458  AliDebug(2, "TClonesArray cleared");
459 
460  fHistStat->Fill(0);
461 
462  // fix for temporary bug in ESDfilter
463  // the AODs with null vertex pointer didn't pass the PhysSel
464  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) return kFALSE;
465 
466  //Event selection
467  Bool_t iseventselected = fCuts->IsEventSelected(fAodEvent);
468  if (!iseventselected) return kFALSE;
469  fHistStat->Fill(1);
470 
471  AliDebug(2, "Event selected");
472 
473  const Int_t nD = fArrayDStartoD0pi->GetEntriesFast();
474  AliDebug(2, Form("Found %d vertices", nD));
475  if (!fUseMCInfo) fHistStat->Fill(2, nD);
476 
477  Int_t pdgMeson = 413;
478  if (fCandidateType == kD0toKpi) pdgMeson = 421;
479 
480  fNCand = 0;
481  fNSBCand = 0;
482 
483  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
484  Int_t isSelected = 0;
485 
486  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fArrayDStartoD0pi->At(icharm)); // D candidates
487  if (!charmCand) continue;
488 
489  Int_t nprongs = charmCand->GetNProngs();
490  AliDebug(2, Form("Candidate is %d, and nprongs = %d", fCandidateType, nprongs));
491 
492  AliAODRecoCascadeHF* dstar = 0;
493 
494  if (fCandidateType == kDstartoKpipi) {
495  dstar = dynamic_cast<AliAODRecoCascadeHF*>(charmCand);
496  if (!dstar) {
497  Error("AliAnalysisTaskSEDmesonsFilterCJ::UserExec","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
498  continue;
499  }
500  }
501 
502  // region of interest + cuts
503  if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(pdgMeson))) continue;
504 
506  FillDstarSideBands(dstar);
507  }
508 
509  //candidate selected by cuts and PID
510  isSelected = fCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
511 
512  if (!isSelected) continue;
513 
514  if (fCandidateType == kDstartoKpipi) {
515  ProcessDstar(dstar, isSelected);
516  }
517  else {
518  ProcessD0(charmCand, isSelected);
519  }
520  } // end of D cand loop
521 
522  AliDebug(2, "Loop done");
523 
524  if (fCombineDmesons) {
525  if (fCombinedDmesons->GetEntriesFast() > 0) {
527  }
528 
529  if (fCombinedDmesonsBkg->GetEntriesFast() > 0) {
531  }
532  }
533 
534  fHistNCandEv->Fill(fCandidateArray->GetEntriesFast());
536  Int_t nsbcand = fSideBandArray->GetEntriesFast();
537  fHistStat->Fill(4, nsbcand);
538  fHistNSBCandEv->Fill(nsbcand);
539  }
540 
541  PostData(1, fOutput);
542  PostData(3, fCandidateArray);
543  PostData(4, fSideBandArray);
544  PostData(5, fCombinedDmesons);
545  PostData(6, fCombinedDmesonsBkg);
546 
547  AliDebug(2, "Exiting method");
548 
549  return kTRUE;
550 }
551 
552 //_______________________________________________________________________________
554 {
555  // Process the D0 candidate.
556 
557  // For MC analysis look for the AOD MC particle associated with the charm candidate
558  AliAODMCParticle* charmPart = 0;
559  if (fUseMCInfo) {
560  Int_t mcLabel = charmCand->MatchToMC(421, fMCarray, fNProngs, fPDGdaughters);
561 
562  if (mcLabel >= 0) {
563  charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(mcLabel));
564  }
565  }
566 
567  if (charmPart) {
568 
569  Int_t origin = CheckOrigin(charmPart, fMCarray);
570  if (origin < 0) return;
571 
572  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
573  fHistStat->Fill(6);
574  return;
575  }
576  if (fRejectDfromB && origin == kFromBottom) {
577  fHistStat->Fill(7);
578  return;
579  }
580  if (fKeepOnlyDfromB && origin != kFromBottom) {
581  fHistStat->Fill(8);
582  return;
583  }
584 
585  fHistStat->Fill(2);
586  }
587  else {
588  fHistStat->Fill(5);
589  }
590 
591  // For MC background fill fSideBandArray
592  if (fUseMCInfo && !charmPart) {
593  if (fUseReco) {
594  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoDecayHF2Prong(*charmCand);
595  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
596 
597  if (fCombineDmesons) {
598  if(!fMultCand) new ((*fCombinedDmesonsBkg)[fNSBCand]) AliEmcalParticle(charmCand);
599  else if(fNSBCand==fAnalyseCand) new ((*fCombinedDmesonsBkg)[0]) AliEmcalParticle(charmCand);
600  }
601 
602  fHistImpParB->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
603  fHistImpParB->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
604  fHistInvMassB->Fill(charmCand->InvMassD0());
605  fHistInvMassB->Fill(charmCand->InvMassD0bar());
606 
607  fNSBCand++;
608  }
609  }
610  // For data or MC signal fill fCandidateArray
611  else {
612  // For data or MC with the requirement fUseReco fill with candidates
613  if (fUseReco) {
614  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODRecoDecayHF2Prong(*charmCand);
615  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
616 
617  if (fCombineDmesons) {
618  if(!fMultCand) new ((*fCombinedDmesons)[fNCand]) AliEmcalParticle(charmCand);
619  else if(fNCand==fAnalyseCand) new ((*fCombinedDmesons)[0]) AliEmcalParticle(charmCand);
620  }
621 
622  fHistImpParS->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
623  fHistImpParS->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
624  fHistInvMassS->Fill(charmCand->InvMassD0());
625  fHistInvMassS->Fill(charmCand->InvMassD0bar());
626  }
627  // For MC with requirement particle level fill with AliAODMCParticle
628  else {
629  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODMCParticle(*charmPart);
630  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODMCParticle(*charmPart);
631  }
632  fHistStat->Fill(3);
633  fNCand++;
634  }
635 
636  // Now filling some more histograms
637 
638  // mass vs pt
639  if (isSelected == 1 || isSelected == 3) fHistInvMassPtD->Fill(charmCand->InvMassD0(), charmCand->Pt());
640  if (isSelected >= 2) fHistInvMassPtD->Fill(charmCand->InvMassD0bar(), charmCand->Pt());
641 
642  if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
643  Int_t isD0 = 0;
644  if (charmPart) {
645  Int_t pdgCode = charmPart->GetPdgCode();
646 
647  if (pdgCode == 421) {
648  isD0 = 1;
649  }
650  else if (pdgCode == -421) {
651  isD0 = -1;
652  }
653  else {
654  AliDebug(2, "Not a D0/D0bar meson!");
655  return;
656  }
657  }
658 
659  FillD0MCTruthKinHistos(charmCand, isSelected, isD0);
660  }
661 }
662 
663 //_______________________________________________________________________________
665 {
666  // Process the D* candidate.
667 
668  AliDebug(2, "Entering method");
669 
670  // For MC analysis look for the AOD MC particle associated with the charm candidate
671  AliAODMCParticle* charmPart = 0;
672  if (fUseMCInfo) {
673  //D* and D0 prongs needed to MatchToMC method
674  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
675  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
676 
677  Int_t mcLabel = dstar->MatchToMC(413, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCarray);
678 
679  if (mcLabel >= 0) {
680  charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(mcLabel));
681  }
682  }
683 
684  if (charmPart) {
685 
686  Int_t origin = CheckOrigin(charmPart, fMCarray);
687  if (origin < 0) return;
688 
689  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
690  fHistStat->Fill(6);
691  return;
692  }
693  if (fRejectDfromB && origin == kFromBottom) {
694  fHistStat->Fill(7);
695  return;
696  }
697  if (fKeepOnlyDfromB && origin != kFromBottom) {
698  fHistStat->Fill(8);
699  return;
700  }
701 
702  fHistStat->Fill(2);
703  }
704  else {
705  fHistStat->Fill(5);
706  }
707 
708  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
709 
710  // For MC background fill fSideBandArray
711  if (fUseMCInfo && !charmPart) {
712  if (fUseReco) {
713  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoCascadeHF(*dstar);
714  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoCascadeHF(*dstar);
715 
716  if (fCombineDmesons) {
717  if(!fMultCand) new ((*fCombinedDmesonsBkg)[fNSBCand]) AliEmcalParticle(dstar);
718  else if(fNSBCand==fAnalyseCand) new ((*fCombinedDmesonsBkg)[0]) AliEmcalParticle(dstar);
719  }
720 
721  fHistInvMassB->Fill(dstar->DeltaInvMass());
722  fHistImpParB->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
723  fHistImpParB->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
724  fHistImpParB->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
725 
726  fNSBCand++;
727  }
728  }
729  // For data and MC signal fill fCandidateArray
730  else {
731  // For data or MC signal with the requirement fUseReco fill with candidates
732  if (fUseReco) {
733  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODRecoCascadeHF(*dstar);
734  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODRecoCascadeHF(*dstar);
735 
736  if (fCombineDmesons) {
737  if(!fMultCand) new ((*fCombinedDmesons)[fNCand]) AliEmcalParticle(dstar);
738  else if(fNCand==fAnalyseCand) new ((*fCombinedDmesons)[0]) AliEmcalParticle(dstar);
739  }
740  }
741  // For MC signal with requirement particle level fill with AliAODMCParticle
742  else {
743  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODMCParticle(*charmPart);
744  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODMCParticle(*charmPart);
745  }
746  fNCand++;
747 
748  fHistInvMassS->Fill(dstar->DeltaInvMass());
749  fHistImpParS->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
750  fHistImpParS->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
751  fHistImpParS->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
752  fHistStat->Fill(3);
753  }
754 
755  // Now filling some more histograms.
756 
757  // select D* in the D0 window.
758  // In the cut object window is loose to allow for side bands
759 
760  // retrieve the corresponding pt bin for the candidate
761  // and set the expected D0 width (x3)
762 
763  Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
764  Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
765  Double_t invMassPDG = mPDGDstar - mPDGD0;
766 
767  Double_t ptD = dstar->Pt();
768  Int_t ptDbin = fCuts->PtBin(ptD);
769  if (ptDbin < 0 || ptDbin >= fCuts->GetNPtBins()) {
770  AliError(Form("Pt %.3f out of bounds", ptD));
771  return;
772  }
773 
774  AliDebug(1, Form("Pt bin %d and sigma D0 %.4f", ptDbin, fSigmaD0[ptDbin]));
775  //consider the D* candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
776  if ((dstar->InvMassD0()>=(mPDGD0-3.*fSigmaD0[ptDbin])) && (dstar->InvMassD0()<=(mPDGD0+3.*fSigmaD0[ptDbin]))) {
777  AliDebug(2, "D0 mass within 3 sigma");
778 
779  // D* delta mass
780  AliDebug(2, "Filling invariant mass vs pt histogram");
781  fHistInvMassPtD->Fill(dstar->DeltaInvMass(), ptD); // 2 D slice for pt bins
782 
783  // Soft pion pt for good candidates
784  Double_t invmassDelta = dstar->DeltaInvMass();
785  if (TMath::Abs(invmassDelta - invMassPDG) < 0.0021) {
786  AliDebug(2, "Filling pion pt histogram");
787  //softpion from D* decay
788  AliAODTrack *softPionTrack = static_cast<AliAODTrack*>(dstar->GetBachelor());
789  if (softPionTrack) fHistPtPion->Fill(softPionTrack->Pt());
790  }
791  }
792 
793  if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
794  Int_t isDstar = charmPart == NULL ? 0 : 1;
795  FillDStarMCTruthKinHistos(dstar, isDstar, isSelected);
796  }
797 
798  AliDebug(2, "Exiting method");
799 }
800 
801 //_______________________________________________________________________________
803 {
804  // Fill some histogram with kinematic information of the D0 candidate.
805 
806  Double_t ptD = charmCand->Pt();
807  Double_t aD = charmCand->Phi();
808  Double_t adaugh[2] = {charmCand->PhiProng(0), charmCand->PhiProng(1)};
809  AliAODTrack* p0 = static_cast<AliAODTrack*>(charmCand->GetDaughter(0));
810  AliAODTrack* p1 = static_cast<AliAODTrack*>(charmCand->GetDaughter(1));
811  Float_t dR0 = DeltaR(charmCand, p0);
812  Float_t dR1 = DeltaR(charmCand, p1);
813 
814  if (isD0 == 0) { //background
815  if (isSelected == 1 || isSelected == 3) { // selected as D0
816  fHistAlphaDKB->Fill(aD-adaugh[0],ptD);
817  fHistAlphaDpiB->Fill(aD-adaugh[1],ptD);
818 
819  fHistDeltaRDKB->Fill(dR0,ptD);
820  fHistDeltaRDpiB->Fill(dR1,ptD);
821  }
822  if (isSelected >= 2) { //selected as D0bar
823  fHistAlphaDpiB->Fill(aD-adaugh[0],ptD);
824  fHistAlphaDKB->Fill(aD-adaugh[1],ptD);
825 
826  fHistDeltaRDpiB->Fill(dR0,ptD);
827  fHistDeltaRDKB->Fill(dR1,ptD);
828  }
829  }
830  else if (isD0 == 1) { //D0
831  fHistAlphaDKS->Fill(aD-adaugh[0],ptD);
832  fHistAlphaDpiS->Fill(aD-adaugh[1],ptD);
833 
834  fHistDeltaRDKS->Fill(dR0,ptD);
835  fHistDeltaRDpiS->Fill(dR1,ptD);
836 
837  if (isSelected == 3) { // selected as both D0bar/D0
838  fHistAlphaDpiR->Fill(aD-adaugh[0],ptD);
839  fHistAlphaDKR->Fill(aD-adaugh[1],ptD);
840 
841  fHistDeltaRDpiR->Fill(dR0,ptD);
842  fHistDeltaRDKR->Fill(dR1,ptD);
843  }
844  }
845  else if (isD0 == -1) { //D0bar
846  fHistAlphaDKS->Fill(aD-adaugh[1],ptD);
847  fHistAlphaDpiS->Fill(aD-adaugh[0],ptD);
848 
849  fHistDeltaRDKS->Fill(dR1,ptD);
850  fHistDeltaRDpiS->Fill(dR0,ptD);
851 
852  if (isSelected == 3) { // selected as both D0bar/D0
853  fHistAlphaDpiR->Fill(aD-adaugh[1],ptD);
854  fHistAlphaDKR->Fill(aD-adaugh[0],ptD);
855 
856  fHistDeltaRDpiR->Fill(dR1, ptD);
857  fHistDeltaRDKR->Fill(dR0, ptD);
858  }
859  }
860 }
861 
862 //_______________________________________________________________________________
864 {
865  // Fill some histogram with kinematic information of the D0 candidate.
866 
867  AliAODTrack *softPionTrack = static_cast<AliAODTrack*>(dstar->GetBachelor());
868  Double_t ptD = dstar->Pt();
869  Double_t aD = dstar->Phi();
870  Double_t apis= softPionTrack->Phi();
871 
872  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
873  Double_t aD0 = D0fromDstar->Phi();
874  Int_t isD0 = D0fromDstar->Charge()>0 ? kTRUE : kFALSE;
875  Double_t aK = isD0 ? D0fromDstar->PhiProng(0) : D0fromDstar->PhiProng(1);
876  Double_t api = isD0 ? D0fromDstar->PhiProng(1) : D0fromDstar->PhiProng(0);
877  Double_t dRDD0 = DeltaR(dstar,D0fromDstar);
878  Double_t dRDpis = DeltaR(dstar,softPionTrack);
879  Double_t dRDpi = DeltaR(dstar, isD0 ? static_cast<AliVParticle*>(D0fromDstar->GetDaughter(1)) : static_cast<AliVParticle*>(D0fromDstar->GetDaughter(0)));
880  Double_t dRDK = DeltaR(dstar, isD0 ? static_cast<AliVParticle*>(D0fromDstar->GetDaughter(0)) : static_cast<AliVParticle*>(D0fromDstar->GetDaughter(1)));
881 
882  if (isDstar) {
883  fHistAlphaDDS ->Fill(aD-aD0, ptD);
884  fHistAlphaDpisS->Fill(aD-apis, ptD);
885  fHistAlphaDpiS ->Fill(aD-api, ptD);
886  fHistAlphaDKS ->Fill(aD-aK, ptD);
887 
888  fHistDeltaRDDS ->Fill(dRDD0, ptD);
889  fHistDeltaRDpisS->Fill(dRDpis, ptD);
890  fHistDeltaRDpiS ->Fill(dRDpi, ptD);
891  fHistDeltaRDKS ->Fill(dRDK, ptD);
892  }
893  else {
894  fHistAlphaDDB ->Fill(aD-aD0, ptD);
895  fHistAlphaDpisB->Fill(aD-apis, ptD);
896  fHistAlphaDpiB ->Fill(aD-api, ptD);
897  fHistAlphaDKB ->Fill(aD-aK, ptD);
898 
899  fHistDeltaRDDB ->Fill(dRDD0, ptD);
900  fHistDeltaRDpisB->Fill(dRDpis, ptD);
901  fHistDeltaRDpiB ->Fill(dRDpi, ptD);
902  fHistDeltaRDKB ->Fill(dRDK, ptD);
903  }
904 }
905 
906 //_______________________________________________________________________________
908 {
909  // Fills the array of Dstar side band candidates.
910 
911  //select by track cuts the side band candidates (don't want mass cut)
912  Int_t isSelected = fCuts->IsSelected(dstar, AliRDHFCuts::kTracks, fAodEvent);
913  if (!isSelected) return;
914 
915  //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
916  Int_t bin = fCuts->PtBin(dstar->Pt());
917  if (bin < 0 || bin >= fCuts->GetNPtBins()) {
918  AliError(Form("Pt %.3f out of bounds", dstar->Pt()));
919  return;
920  }
921 
922  const Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
923 
924  //if data and Dstar from D0 side band
925  if (((dstar->InvMassD0() < (mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0() > (mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/ ||
926  ((dstar->InvMassD0() > (mPDGD0+3.*fSigmaD0[bin])) && (dstar->InvMassD0() < (mPDGD0+10.*fSigmaD0[bin]))) /*right side band*/) {
927 
928  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoCascadeHF(*dstar);
929  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoCascadeHF(*dstar);
930  fNSBCand++;
931  }
932 }
933 
934 //_______________________________________________________________________________
936 {
937  // Set the mass limits using the PDG code and the given range.
938 
939  Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
940 
941  // compute the Delta mass for the D*
942  if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
943 
944  fMinMass = mass - range;
945  fMaxMass = mass + range;
946 
947  AliInfo(Form("Setting mass limits to %f, %f", fMinMass, fMaxMass));
948  if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) {
949  AliError(Form("Wrong mass limits! Task '%s' will be disabled!", GetName()));
950  fInhibitTask = kTRUE;
951  }
952 }
953 
954 //_______________________________________________________________________________
956 {
957  // Set the mass limits.
958 
959  if (uplimit>lowlimit) {
960  fMinMass = lowlimit;
961  fMaxMass = uplimit;
962  }
963  else {
964  printf("Error! Lower limit larger than upper limit!\n");
965  fMinMass = uplimit - uplimit*0.2;
966  fMaxMass = uplimit;
967  }
968 }
969 
970 //_______________________________________________________________________________
972 {
973  // Set the D0 width for the D* analysis.
974 
975  if (nptbins > 30) {
976  AliWarning("Maximum number of bins allowed is 30!");
977  return kFALSE;
978  }
979 
980  if (!width) return kFALSE;
981  for (Int_t ipt=0; ipt < nptbins; ipt++) fSigmaD0[ipt] = width[ipt];
982 
983  return kTRUE;
984 }
985 
986 //_______________________________________________________________________________
988 {
989  // Allocate the histograms for the analysis.
990 
991  // Statistics
992  fHistStat = new TH1I("fHistStat", "Statistics", 9, -0.5, 8.5);
993  fHistStat->GetXaxis()->SetBinLabel(1, "N ev anal");
994  fHistStat->GetXaxis()->SetBinLabel(2, "N ev sel");
995  if (fUseMCInfo) {
996  if (fUseReco) {
997  fHistStat->GetXaxis()->SetBinLabel(3, "N D");
998  }
999  else {
1000  fHistStat->GetXaxis()->SetBinLabel(3, "N Gen D");
1001  }
1002  }
1003  else {
1004  fHistStat->GetXaxis()->SetBinLabel(3, "N Cand");
1005  }
1006  fHistStat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
1007  if (fCandidateType == kDstartoKpipi) {
1008  fHistStat->GetXaxis()->SetBinLabel(5, "N side band cand");
1009  }
1010  if (fUseMCInfo) {
1011  fHistStat->GetXaxis()->SetBinLabel(6, "N Background");
1012  fHistStat->GetXaxis()->SetBinLabel(7, "N rej no quark");
1013  fHistStat->GetXaxis()->SetBinLabel(8, "N rej from B");
1014  fHistStat->GetXaxis()->SetBinLabel(9, "N rej from D");
1015  }
1016 
1017  fHistStat->SetNdivisions(1);
1018  fOutput->Add(fHistStat);
1019 
1020  fHistNCandEv = new TH1F("fHistNCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0., 100.);
1021  fOutput->Add(fHistNCandEv);
1022 
1023  if ((fCandidateType == kDstartoKpipi) || fUseMCInfo) {
1024  fHistNSBCandEv = new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
1025  fOutput->Add(fHistNSBCandEv);
1026  }
1027 
1028  if (fCandidateType == kDstartoKpipi) {
1029  fHistPtPion = new TH1F("fHistPtPion", "Primary pions candidates pt", 500, 0., 10.);
1030  fHistPtPion->SetStats(kTRUE);
1031  fHistPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1032  fHistPtPion->GetYaxis()->SetTitle("entries");
1033  fOutput->Add(fHistPtPion);
1034  }
1035 
1036  // Invariant mass related histograms
1037  const Int_t nbinsmass = 200;
1038  const Int_t ptbinsD = 100;
1039  Float_t ptmin =0.;
1040  Float_t ptmax =50.;
1041  fHistInvMassPtD = new TH2F("fHistInvMassPtD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, ptbinsD, ptmin, ptmax);
1042  fHistInvMassPtD->SetStats(kTRUE);
1043  fHistInvMassPtD->GetXaxis()->SetTitle("mass (GeV/c)");
1044  fHistInvMassPtD->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1045  fOutput->Add(fHistInvMassPtD);
1046 
1047  fHistImpParS = new TH2F("fHistImpParS", "Impact parameter of daughter tracks; Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
1048  fOutput->Add(fHistImpParS);
1049 
1050  fHistInvMassS = new TH1F("fHistInvMassS", "D invariant mass distribution (filled with fCandidateArray)", nbinsmass, fMinMass, fMaxMass);
1051  fHistInvMassS->SetStats(kTRUE);
1052  fHistInvMassS->GetXaxis()->SetTitle("mass (GeV/c)");
1053  fHistInvMassS->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1054  fOutput->Add(fHistInvMassS);
1055 
1056  const Int_t nbinsalpha = 200;
1057  Float_t minalpha = -TMath::Pi();
1058  Float_t maxalpha = TMath::Pi();
1059  const Int_t nbinsdeltaR = 200;
1060  Float_t mindeltaR = 0.;
1061  Float_t maxdeltaR = 10.;
1062 
1063  if (fUseMCInfo) {
1064  fHistImpParB = new TH2F("fHistImpParB", "Impact parameter of daughter tracks (Background); Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
1065  fOutput->Add(fHistImpParB);
1066 
1067  fHistInvMassB = new TH1F("fHistInvMassB", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass);
1068  fHistInvMassB->SetStats(kTRUE);
1069  fHistInvMassB->GetXaxis()->SetTitle("mass (GeV/c)");
1070  fHistInvMassB->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1071  fOutput->Add(fHistInvMassB);
1072 
1073  if (fCandidateType == kDstartoKpipi) {
1074 
1075  fHistAlphaDDS = new TH2F("fHistAlphaDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1076  fHistAlphaDpisS = new TH2F("fHistAlphaDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1077  fHistAlphaDpiS = new TH2F("fHistAlphaDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1078  fHistAlphaDKS = new TH2F("fHistAlphaDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1079 
1080  fHistAlphaDDB = new TH2F("fHistAlphaDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1081  fHistAlphaDpisB = new TH2F("fHistAlphaDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1082  fHistAlphaDpiB = new TH2F("fHistAlphaDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1083  fHistAlphaDKB = new TH2F("fHistAlphaDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1084 
1085  fHistDeltaRDDS = new TH2F("fHistDeltaRDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1086  fHistDeltaRDpisS = new TH2F("fHistDeltaRDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1087  fHistDeltaRDpiS = new TH2F("fHistDeltaRDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1088  fHistDeltaRDKS = new TH2F("fHistDeltaRDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1089 
1090  fHistDeltaRDDB = new TH2F("fHistDeltaRDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1091  fHistDeltaRDpisB = new TH2F("fHistDeltaRDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1092  fHistDeltaRDpiB = new TH2F("fHistDeltaRDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1093  fHistDeltaRDKB = new TH2F("fHistDeltaRDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1094 
1095  fOutput->Add(fHistAlphaDDS);
1096  fOutput->Add(fHistAlphaDpisS);
1097  fOutput->Add(fHistAlphaDpiS);
1098  fOutput->Add(fHistAlphaDKS);
1099  fOutput->Add(fHistAlphaDDB);
1100  fOutput->Add(fHistAlphaDpisB);
1101  fOutput->Add(fHistAlphaDpiB);
1102  fOutput->Add(fHistAlphaDKB);
1103 
1104  fOutput->Add(fHistDeltaRDDS);
1105  fOutput->Add(fHistDeltaRDpisS);
1106  fOutput->Add(fHistDeltaRDpiS);
1107  fOutput->Add(fHistDeltaRDKS);
1108  fOutput->Add(fHistDeltaRDDB);
1109  fOutput->Add(fHistDeltaRDpisB);
1110  fOutput->Add(fHistDeltaRDpiB);
1111  fOutput->Add(fHistDeltaRDKB);
1112  }
1113  else if (fCandidateType == kD0toKpi) {
1114 
1115  fHistAlphaDpiS = new TH2F("fHistAlphaDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1116  fHistAlphaDKS = new TH2F("fHistAlphaDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1117  fHistAlphaDpiR = new TH2F("fHistAlphaDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1118  fHistAlphaDKR = new TH2F("fHistAlphaDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1119 
1120  fHistAlphaDpiB = new TH2F("fHistAlphaDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1121  fHistAlphaDKB = new TH2F("fHistAlphaDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1122 
1123  fHistDeltaRDpiS = new TH2F("fHistDeltaRDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1124  fHistDeltaRDKS = new TH2F("fHistDeltaRDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1125  fHistDeltaRDpiR = new TH2F("fHistDeltaRDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1126  fHistDeltaRDKR = new TH2F("fHistDeltaRDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1127 
1128  fHistDeltaRDpiB = new TH2F("fHistDeltaRDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1129  fHistDeltaRDKB = new TH2F("fHistDeltaRDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1130 
1131  fOutput->Add(fHistAlphaDpiS);
1132  fOutput->Add(fHistAlphaDKS);
1133  fOutput->Add(fHistAlphaDpiR);
1134  fOutput->Add(fHistAlphaDKR);
1135  fOutput->Add(fHistAlphaDpiB);
1136  fOutput->Add(fHistAlphaDKB);
1137 
1138  fOutput->Add(fHistDeltaRDpiS);
1139  fOutput->Add(fHistDeltaRDKS);
1140  fOutput->Add(fHistDeltaRDpiR);
1141  fOutput->Add(fHistDeltaRDKR);
1142  fOutput->Add(fHistDeltaRDpiB);
1143  fOutput->Add(fHistDeltaRDKB);
1144  }
1145  }
1146 
1147  return kTRUE;
1148 }
1149 
1150 //_______________________________________________________________________________
1151 Float_t AliAnalysisTaskSEDmesonsFilterCJ::DeltaR(AliVParticle *p1, AliVParticle *p2) const
1152 {
1153  // Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1154 
1155  if (!p1 || !p2) return -1;
1156 
1157  Double_t phi1 = p1->Phi();
1158  Double_t eta1 = p1->Eta();
1159  Double_t phi2 = p2->Phi();
1160  Double_t eta2 = p2->Eta();
1161 
1162  Double_t dPhi = phi1 - phi2;
1163  if(dPhi < -(TMath::Pi())/2) dPhi = dPhi + TMath::TwoPi();
1164  if(dPhi > (3*(TMath::Pi()))/2) dPhi = dPhi - TMath::TwoPi();
1165 
1166  Double_t dEta = eta1 - eta2;
1167 
1168  Double_t deltaR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
1169 
1170  return deltaR;
1171 }
1172 
1173 //_______________________________________________________________________________
1175 {
1176  //
1177  // Add event tracks to a collection that already contains the D candidates, excluding the daughters of the D candidates
1178  //
1179 
1180  if (!tracks) return;
1181 
1182  TObjArray allDaughters(10);
1183  allDaughters.SetOwner(kFALSE);
1184 
1185  TIter next(coll);
1186  AliEmcalParticle* emcpart = 0;
1187  while ((emcpart = static_cast<AliEmcalParticle*>(next()))) {
1188  AliAODRecoDecay* reco = dynamic_cast<AliAODRecoDecay*>(emcpart->GetTrack());
1189  AliDebug(2, Form("Found a D meson candidtate with pT = %.3f, eta = %.3f, phi = %.3f", reco->Pt(), reco->Eta(), reco->Phi()));
1190  if (reco) AddDaughters(reco, allDaughters);
1191  }
1192 
1193  tracks->ResetCurrentID();
1194  AliVTrack* track = 0;
1195  Int_t n = coll->GetEntriesFast();
1196  while ((track = static_cast<AliVTrack*>(tracks->GetNextAcceptParticle()))) {
1197  if (allDaughters.Remove(track) == 0) {
1198  new ((*coll)[n]) AliEmcalParticle(track);
1199  n++;
1200  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is included", tracks->GetCurrentID(), track->Pt(), track->Eta(), track->Phi()));
1201  }
1202  else {
1203  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is excluded", tracks->GetCurrentID(), track->Pt(), track->Eta(), track->Phi()));
1204  }
1205  }
1206 }
1207 
1208 //_______________________________________________________________________________
1210 {
1211  // Add all the dauthers of cand in an array. Follows all the decay cascades.
1212 
1213  Int_t n = cand->GetNDaughters();
1214 
1215  //Printf("AddDaughters: the number of dauhters is %d", n);
1216 
1217  Int_t ntot = 0;
1218  Double_t pt = 0;
1219  for (Int_t i = 0; i < n; i++) {
1220  AliVTrack* track = dynamic_cast<AliVTrack*>(cand->GetDaughter(i));
1221  if (!track) continue;
1222 
1223  AliAODRecoDecay* cand2 = dynamic_cast<AliAODRecoDecay*>(track);
1224 
1225  if (cand2) {
1226  //Printf("Daughter pT = %.3f --> ", track->Pt());
1227  pt += AddDaughters(cand2, daughters);
1228  }
1229  else {
1230  if (!track->InheritsFrom("AliAODTrack")) {
1231  Printf("Warning: One of the daughters is not of type 'AliAODTrack' nor 'AliAODRecoDecay'.");
1232  continue;
1233  }
1234  //Printf("Daughter pT = %.3f", track->Pt());
1235  daughters.AddLast(track);
1236  pt += track->Pt();
1237  ntot++;
1238  }
1239  }
1240 
1241  //Printf("Total pt of the daughters = %.3f", pt);
1242 
1243  return pt;
1244 }
1245 
1246 //_________________________________________________________________________________________________
1247 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODRecoDecay* cand, TClonesArray* mcArray)
1248 {
1249  // Checks whether the mother of the D meson candidate comes from a charm or a bottom quark.
1250 
1251  if (!mcArray) return -1;
1252 
1253  Int_t labDau0 = static_cast<AliVTrack*>(cand->GetDaughter(0))->GetLabel();
1254  if (labDau0 < 0) return -1;
1255 
1256  AliAODMCParticle* part = static_cast<AliAODMCParticle*>(mcArray->At(labDau0));
1257  return CheckOrigin(part, mcArray);
1258 }
1259 
1260 //_________________________________________________________________________________________________
1261 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODRecoDecay* cand, AliStack* stack)
1262 {
1263  // Checks whether the mother of the D meson candidate comes from a charm or a bottom quark.
1264 
1265  if (!stack) return -1;
1266 
1267  Int_t labDau0 = static_cast<AliVTrack*>(cand->GetDaughter(0))->GetLabel();
1268  if (labDau0 < 0) return -1;
1269 
1270  return CheckOrigin(labDau0, stack);
1271 }
1272 
1273 //_________________________________________________________________________________________________
1274 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODMCParticle* part, TClonesArray* mcArray)
1275 {
1276  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1277 
1278  if (!part) return -1;
1279  if (!mcArray) return -1;
1280 
1281  Int_t pdgGranma = 0;
1282  Int_t mother = part->GetMother();
1283  Int_t istep = 0;
1284  Int_t abspdgGranma = 0;
1285  Bool_t isFromB = kFALSE;
1286  Bool_t isQuarkFound = kFALSE;
1287 
1288  while (mother >= 0) {
1289  istep++;
1290  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1291  if (mcGranma >= 0) {
1292  pdgGranma = mcGranma->GetPdgCode();
1293  abspdgGranma = TMath::Abs(pdgGranma);
1294  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1295  isFromB = kTRUE;
1296  }
1297 
1298  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1299  mother = mcGranma->GetMother();
1300  }
1301  else {
1302  ::Error("AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1303  break;
1304  }
1305  }
1306 
1307  if (isQuarkFound) {
1308  if (isFromB) {
1309  return kFromBottom;
1310  }
1311  else {
1312  return kFromCharm;
1313  }
1314  }
1315  else {
1316  return kQuarkNotFound;
1317  }
1318 }
1319 
1320 //_________________________________________________________________________________________________
1322 {
1323  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1324 
1325  if (!stack) return -1;
1326 
1327  TParticle* part = stack->Particle(ipart);
1328  if (!part) return -1;
1329 
1330  Int_t pdgGranma = 0;
1331  Int_t mother = part->GetFirstMother();
1332  Int_t istep = 0;
1333  Int_t abspdgGranma = 0;
1334  Bool_t isFromB = kFALSE;
1335  Bool_t isQuarkFound = kFALSE;
1336 
1337  while (mother >= 0) {
1338  istep++;
1339  TParticle* mcGranma = stack->Particle(mother);
1340  if (mcGranma >= 0) {
1341  pdgGranma = mcGranma->GetPdgCode();
1342  abspdgGranma = TMath::Abs(pdgGranma);
1343  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1344  isFromB = kTRUE;
1345  }
1346 
1347  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1348  mother = mcGranma->GetFirstMother();
1349  }
1350  else {
1351  ::Error("AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1352  break;
1353  }
1354  }
1355 
1356  if (isQuarkFound) {
1357  if (isFromB) {
1358  return kFromBottom;
1359  }
1360  else {
1361  return kFromCharm;
1362  }
1363  }
1364  else {
1365  return kQuarkNotFound;
1366  }
1367 }
1368 
1369 //_________________________________________________________________________________________________
1370 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckDecayChannel(AliAODMCParticle* part, TClonesArray* mcArray)
1371 {
1372  // Determine the decay channel
1373 
1374  if (!part) return -1;
1375  if (!mcArray) return -1;
1376 
1378 
1379  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1380 
1381  if (part->GetNDaughters() == 2) {
1382 
1383  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1384  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1385 
1386  if (!d1 || !d2) {
1387  return decay;
1388  }
1389 
1390  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1391  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1392 
1393  if (absPdgPart == 421) { // D0 -> K pi
1394 
1395  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1396  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1397  decay = kDecayD0toKpi;
1398  }
1399  }
1400 
1401  if (absPdgPart == 413) { // D* -> D0 pi
1402 
1403  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1404  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1405  if (D0decay == kDecayD0toKpi) {
1406  decay = kDecayDStartoKpipi;
1407  }
1408  }
1409 
1410  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1411  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1412  if (D0decay == kDecayD0toKpi) {
1413  decay = kDecayDStartoKpipi;
1414  }
1415  }
1416  }
1417  }
1418 
1419  return decay;
1420 }
1421 
1422 //_________________________________________________________________________________________________
1424 {
1425  // Determine the decay channel
1426 
1427  if (!stack) return -1;
1428 
1429  TParticle* part = stack->Particle(ipart);
1430 
1431  if (!part) return -1;
1432 
1434 
1435  if (part->GetNDaughters() == 2) {
1436 
1437  Int_t id1 = part->GetDaughter(0);
1438  Int_t id2 = part->GetDaughter(1);
1439 
1440  TParticle* d1 = stack->Particle(id1);
1441  TParticle* d2 = stack->Particle(id2);
1442 
1443  if (!d1 || !d2) {
1444  return decay;
1445  }
1446 
1447  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1448  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1449 
1450 
1451  if (part->GetPdgCode() == 421) { // D0 -> K pi
1452 
1453  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1454  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1455  decay = kDecayD0toKpi;
1456  }
1457  }
1458 
1459  if (part->GetPdgCode() == 413) { // D* -> D0 pi
1460 
1461  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1462  Int_t D0decay = CheckDecayChannel(id1, stack);
1463  if (D0decay == kDecayD0toKpi) {
1464  decay = kDecayDStartoKpipi;
1465  }
1466  }
1467 
1468  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1469  Int_t D0decay = CheckDecayChannel(id2, stack);
1470  if (D0decay == kDecayD0toKpi) {
1471  decay = kDecayDStartoKpipi;
1472  }
1473  }
1474  }
1475  }
1476 
1477  return decay;
1478 
1479 }
void ProcessDstar(AliAODRecoCascadeHF *dstar, Int_t isSelected)
Int_t pdg
static Int_t CheckOrigin(AliAODRecoDecay *cand, TClonesArray *mcArray)
virtual AliVParticle * GetNextAcceptParticle()
double Double_t
Definition: External.C:58
Definition: External.C:236
TClonesArray * fCombinedDmesonsBkg
contains candidates selected by AliRDHFCuts and the rest of the event tracks
Double_t DeltaInvMass() const
void ProcessD0(AliAODRecoDecayHF2Prong *charmCand, Int_t isSelected)
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
Base task in the EMCAL framework.
Double_t mass
void FillDStarMCTruthKinHistos(AliAODRecoCascadeHF *dstar, Int_t, Int_t isDstar)
void FillD0MCTruthKinHistos(AliAODRecoDecayHF2Prong *charmCand, Int_t isSelected, Int_t isD0)
static Double_t AddDaughters(AliAODRecoDecay *cand, TObjArray &daughters)
void AddEventTracks(TClonesArray *coll, AliParticleContainer *tracks)
AliVTrack * GetTrack() const
Container for particles within the EMCAL framework.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Int_t fNCand
contains bkg candidates selected by AliRDHFCuts and the rest of the event tracks
void FillDstarSideBands(AliAODRecoCascadeHF *dstar)
TClonesArray * fSideBandArray
contains candidates selected by AliRDHFCuts
int Int_t
Definition: External.C:63
Definition: External.C:204
static Int_t CheckDecayChannel(AliAODMCParticle *part, TClonesArray *mcArray)
float Float_t
Definition: External.C:68
const Double_t ptmax
AliAODTrack * GetBachelor() const
const Double_t ptmin
Bool_t SetD0WidthForDStar(Int_t nptbins, Float_t *width)
ClassImp(AliAnalysisTaskSEDmesonsFilterCJ) AliAnalysisTaskSEDmesonsFilterCJ
Bool_t IsEventSelected(AliVEvent *event)
decay
Definition: HFPtSpectrum.C:41
AliEmcalList * fOutput
!output list
Int_t fNSBCand
number of selected D candidates already added to fCandidateArray
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:274
Double_t InvMassD0() const
TH1 * fHistStat
number of selected side-band D candidates already added to fSideBandArray
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:233
Float_t DeltaR(AliVParticle *p1, AliVParticle *p2) const
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
bool Bool_t
Definition: External.C:53
AliAODRecoDecayHF2Prong * Get2Prong() const
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:294
Int_t PtBin(Double_t pt) const
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Int_t nptbins
TClonesArray * fCombinedDmesons
contains candidates selected by AliRDHFCuts::IsSelected(kTracks), to be used for side bands (DStar ca...